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Abstract 

We study two superconducting systems using the Landau-Ginzburg 
equations. The first is a superconducting half-space with an appUed 
magnetic field parallel to the surface. We calculate the maximum 
applied field that still supports superconductivity in the material by 
solving the Landau-Ginzburg equations analytically in asymptotic 
regimes of the Landau-Ginzburg parameter, k. These results are 
then checked against numerical studies. 

The other system is a thin, current-carrying strip. Here we do the 
first systematic study of nucleation of superconductivity in a nor- 
mal strip and of the motion of a normal-superconducting interface 
in the presence of a current. This dissipative system requires solv- 
ing the time-dependent Landau-Ginzburg equations analytically 
and numerically. 
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Chapter 1 

Introduction — Landau-Ginzburg 
Equations 

We are going to examine superconductors in large applied fields that push them to the verge of nonequilibrium 
or out of equilibrium. We are interested in the nonlinear dynamics at a mesoscopic level. The mechanisms 
of superconductivity are not intrinsic to this study. Of more interest is how we can describe balances among 
populations of superconducting and normal electrons or the magnetic field and the Meissner effect. 

We study two systems with phase boundaries. One is the critical field of a first-order phase transition, 
the other the progression of a phase transition in time. One first analyzes such systems in the bulk by looking 
at free energies and possibilities of metastable states, but studying spatial variations in the applied fields 
and responses of the metallic states yields valuable insight into how systems either maintain equilibrium or 
temper their tendency to collapse. 

The Landau-Ginzburg (GL) equations will be our model for these systems. It is a venerable mean-field 
theory which continues to provide insight into the actual behavior of metals. These equations describe general 
phase-transitions in terms of spatially-dependent variables. The Landau-Ginzburg equations describe very 
nonlinear, and hence complex, systems with little price in terms of calculating specific material parameters.[J 

The Landau-Ginzburg equations, as we work with them, are coupled nonlinear differential equations. 
All of the new research in this dissertation deals directly with those differential equations. The focus of 
the work, therefore, is essentially mathematical and of interest in its own right. All of it, however, is 
introduced, executed, and summarized in light of the physical content. In fact, we begin this introduction 
to the Landau-Ginzburg equations with a description of the free energy. 

1.1 Thermodynamic Free Energy 

It can be shown from thermodynamic arguments that the transition of a metal from normal to supercon- 
ducting is first order if the system is in an applied magnetic field. When the metal is normal, the magnetic 
field penetrates the sample. The Gibbs free energy is the free energy of the metal plus that of the magnetic 
field 

G„ = J^n + f^. (1.1) 

^The recent tendency in literature to use the name Landau-Ginzburg equations rather than the traditional Ginzburg-Landau 
equations seems to indicate an interest in ascribing primacy to Landau. In either case, it has been nearly impossible to avoid 
abbreviating the name with either GL or TDGL (for time-dependent Landau-Ginzburg equations.) 
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When the system becomes superconducting, it expels the magnetic field. The additional work required to 
expel the field is equal to the original field energy 

G,^F, + ^. (1.2) 

Here, Fs is the free energy of the metal in the superconducting state. The latent heat of transition from 
normal to superconducting is thus 

7^ flJ-f 

L = T{Sn-S,)^ H—. (1.3) 

Type-i metals abruptly expel applied magnetic fields at a critical field or temperature so that dH/dT is 
nonzero for all applied fields. This is a first-order phase transition that will exhibit hysteresis. When there 
is no applied field, the transition is second-order. 

We want to look at the mechanics of the first-order phase transition on a mesoscopic scale with the 
Landau-Ginzburg equations. The Landau- Ginzburg free energy permits spatial variations in the free energy 
difference between the normal and superconducting states, G„ — Gs, so that we can find what spatial 
distributions minimize the free energy. 



1.2 Landau-Ginzburg Free Energy 

In this section, we write down the free energy of a superconductor and derive Landau-Ginzburg differential 
equations describing its phases. This derivation is well documented in several works |2^, |7^. Our goal here 
is to demonstrate physical meanings of the nondimensionalized forms used in the rest of the work. The 
Landau-Ginzburg equations are derived from a free energy postulated to be valid near phase transitions 
from a less ordered phase to a more ordered one. While the form of the free energy near a phase transition 
is general, we will examine these equations solely in the context of the transition from superconducting to 
normal. 

We are going to look at a phase transition from a higher temperature, disordered state to a lower 
temperature state with more order and less symmetry. Above an onset temperature called the critical 
temperature, Tc, the system has one available state. Below that temperature, the superconducting state 
with lower free energy becomes accessible as shown in Figure 1.1.^ If we can find how the Free Energy 



of the system depends on its superconducting state, then standard thermodynamic arguments will tell us 
when and how the system is superconducting. We measure how superconducting the system is with an order 
parameter, ip = \tp\e'^^ . When the order parameter is zero, the system is normal. When the order parameter 
is nonzero, the system is superconducting. The density of superconducting electrons is associated with the 
norm of the order parameter, oc IV'P- 

The most significant independent variables are the temperature and order parameter. Because the order 
parameter is small at the onset of the phase transition, we want to write the free energy density as a series 
in the order parameter 

f = fn + aiP + h^l;^ +cij^ +dilj^ + ... (1.4) 

where the total free energy is the free energy of the normal state plus an excess due to superconductivity. 
The system, thus the order parameter, will minimize the free energy. Because the order parameter, tp, is 
complex (isomorphic to the S0(2) symmetry group), the lowest order invariants are ji/'p and IV'I^- No first 

^The Free Energy decreases with increasing temperature as seen from the thermodynamic relation S = — > where A 

is the Helmholtz Free Energy. 
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Figure 1.1: This figure shows the free energy of the normal and superconducting phases as a function 
of temperature. Below the critical temperature, Tc, the system undergoes a second-order phase change 
to the superconducting phase. 



or third-order terms can appear in the free energy. The correct form, to fourth order, is 



(1.5) 



The coefficients a and (3 are constants at a given temperature and depend on the system. There can be a 
free energy minimum for a nonzero order parameter only if a < and /3 > 0. In this case, the free energy 

V^a/P. The 



1.2, and the equihbrium state is ip 



as a function of order parameter is as shown in Figure 
parameters, a and depend on temperature, and they define a maximum value of the order parameter so 
that it varies between ■0 = and ip = v— a//?. Equation 1.5 is the Landau-Ginzburg equation because it 
contains the the basic physics of a bulk system. 

Describing a system with any spatial variation requires modifications to the Landau-Ginzburg free energy 
above to account for the energy cost of bending the order parameter. The lowest order derivative of the 
order parameter gives us 



/=/„+«iv^p + fiV'r+^i?ivv'i'. 



(1.6) 



Lastly, if the metal is in an applied magnetic field, H, we can include the vector potential. A, of the magnetic 
induction, B = V x A, gauge-invariantly to find Ginzburg and Landau's expansion of the free energy density 



1 

2m 



■.^ eA 

'ihV 

c 



(1.7) 



/„ is the free energy of the normal phase. The second to last term is just the energy of the magnetic field 
and last term a familiar representation of kinetic energy, here the kinetic energy of superconducting charge 
carriers. 
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1.3 Time-dependent Landau-Ginzburg Equations 



Because we want to look at a metal at constant temperature, we switch to the Gibbs free energy density. 
The Gibbs free energy is generally used for chemical systems and includes the work required to expel the 
magnetic field from the sample. That energy is of the form J7 = ^ / B • H 



B H 

An 



The total Gibbs free energy for the sample, integrated over all space, is then 

eA 



1 
2m 



-ihV - 



B H 

47r 



(l.E 



(1.9) 



The state variables, (^/i. A), will minimize the free energy in steady state. The minimization conditions can 
be written as functional derivatives 



SG 







and 



SG_ 
SA 



0. 



(1.10) 



Equation 1.10 defines the time-independent Ginzburg-Landau equations if G is the GL free energy. If we 
want to describe nonequilibrium systems, we can construct a relaxational model using the same free energy 



SG_ 

dt S^* 



and 



dA 

'dt 



-T- 



^dG_ 
SA' 



(1.11) 



This form is called the Glauber model or Type A dynamics. The rate of relaxation of the system is controlled 
by r, called an Onsager coefficient. While the order parameter and vector potential are here shown with the 
same relaxation rate, we will later choose separate constants, F, for each. 

The Glauber model does not generate entirely general dynamical equations. For instance, important 
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work on thin superconducting strips uses a model including higher order time derivatives 

1 dip 1 d^ip _ SG 



(1.12) 



The second-order derivative in time permits oscillatory, self-driving solutions. Modern field theories typically 
include a Langevin term describing noise 

? = -r|^ + e (1-13) 

The noise term is necessary in order to include thermally-driven transitions explicitly in the theory. Without 
it, the system would always remain in any metastable state. The pure rclaxational dynamics of the Glauber 
model assume the system will always proceed towards a lower free energy without any self-driving terms or 
noise. 

The plan for the rest of the chapter is first to follow the prescription above to derive a preliminary version 
of the time-dependent Landau-Ginzburg equation then write down the improved and preferred equations 
derived by Gorkov and Eliashberg ||5l|. The point is that, while Gorkov and Eliashberg derived time- 
dependent Landau-Ginzburg equations from microscopies, the results are only a minor modification of those 

derived from the simple relaxational dynamics above. 

SG 



If we compute the two functional derivatives, 4pr and |x (Appendix A.l), we find 



1^ + aV' + m'i' + T^H^"^ - -A)V = (1.14) 
1 ot 2m c 

1 r/_A_ 1 p f p p 

-— + —V X (V X A - H) + - — i;{-ihV + -A)V'* -I- - — V*(«?iV + -A)^j = 0. (1.15) 
r ot An 2mc c 2mc c 

The first equation, describing the motion of the order parameter, is largely described above. We can see in 
this equation that in the bulk where ip does not vary, it will have the value V'o = \/ —a/ (3. In the second 
equation, we can identify first the term V x V x A as the total current generating the local magnetic field. If 
we notice that a stationary system has ^ = 0, we see the last two terms must constitute the supercurrent. 
For instance, one of the boundary conditions on this system is that the total current out of the sample be 
zero, or 

[ihV + -A)tP = (1.16) 

c 

VxA==H. (1.17) 
A more common way to see the second equation written is 

The first term, ^ is related to the total current for a system out of equilibrium. Because it expresses the 
inertia of the superconducting electrons accelerating in the presence of a current, it is called the kinetic 
inductance. 

The equation as it stands describes only supercurrent in a sample. Since we are also interested in samples 
with normal current, we should add an electric field somewhere. If we start from simple E&M, we know that 
the total current is related to the magnetic field by 

An An 

V X H= — J = — (J„+J,). (1.19) 

c c 

Distinguishing the current as either normal or superconducting is called the two-fluid model. It implies 
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there are distinct populations of normal and superconducting electrons. A more current understanding is 
that the norm of the order parameter represents the proportion of electrons participating in Cooper pairing. 
Comparing Eqn. 1.19| with Eqn. LIS, we see the supercurrent is 



ieh 
2m 



-IVI'A. 



(1.20) 



Because the supercurrent does not dissipate energy, the normal current alone is responsible for any elec- 
tric field. The two-fluid model seems to identify T^^dA/dt with the normal current. That term is not 
normal current but the kinetic inductance of the superconducting electrons. Accelerating Cooper pairs to 
an equilibrium superfluid velocity generates an ephemeral voltage in the sample even though there is no 
dissipation. 

If we want to include a normal current in our model, we can look at replacing the kinetic inductance 
with a traditional term of the form J„ = crE. Maxwell's Equations tell us we can always express the electric 
field with a scalar potential (p of the form 



E 



1 dA 

c~dt 



(1.21) 



If we solve that equation for J„ and substitute it into Eqn. 1.18, we find the kinetic inductance appears 
naturally as part of the electric field. Assembling the pieces in the form Jj = J„ -I- Jg, we find 



in 



-V X V X A 



adA 



ihe 
2m 



(^vv^* - V* VV') ■ 



47r 



-V X H 



(1.22) 



In addition, we have found that Maxwell's Equations determine the Onsager Coefficient, F for the relaxation 
equation for the vector potential. While there was a single relaxation constant, F, for the time dependence 
of both the order parameter and magnetic field, we will now set one of them to agree with Maxwell's 
equations. This is just a likely choice, and those who look to Landau-Ginzburg equations for more precise 
correspondence with physical values will often reserve the right to vary independently relaxation coefhcients 
for the order parameter and magnetic field. 

One can find in the above equation a gauge symmetry whereby the equation is unchanged under trans- 
formations of the type 



A ^ A + Vx 



idx 

c dt' 



(1.23) 



In order that both the equation for the magnetic potential and the order parameter be gauge invariant, we 
should add a gauge-invariant term for </) to Eqn. 1.14 to get 



1 

f~dt 



le 
KT 



1 

2m 



{-ihV - -A)V = 0. 



(1.24) 



Our insertion of an electric field by hand has effectively switched to proper canonical variables. We have 
been seeking a sensible form similar to that derived from microscopies by Gorkov and Eliashberg |5l| and 
quoted in Du et al. pl|P| 



+ ie$V - DikV + ieAft/j + (/3|?/;P + aU ^ 

ot 



i^VxVxA = E + 2T 



Itpl^A + - V'VV'*) + V X H 

2e 



(1.25) 
(1.26) 



Without defining the constants, we can say that these equations are, term for term, proportional to those 

^When quoting this equation, I have switched the gauge field, A, to match the sign of my definition, B = V X A. Du et ah 
explained the sign changes by switching the sign of the electric field, a choice not generally welcome to physicists. 
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we have derived. For reasonable choices of constants, they behave hke those derived straight from the 
Landau-Ginzburg free energy. 



The derivations of Eqns. 1.24 and 1.22 may have seemed correct enough with respect to MaxweU's Equa- 
tions that the version of the time-dependent Landau-Ginzburg equations derived by Gorkov and Ehashberg 
have an unnecessary number of parameters. A quick way to see how we have been duped by treading 
well-worn paths is to look for a continuity equation for the charge associated with the superconducting 
current 

^ + V.J. = 0. (1.27) 
at 

A note by Neu |8^ points out that such a continuity equation exists, but, in our variables, we would need 
to choose r = —ie/h. That choice of the Onsager coefficient would violate our relaxational dynamics, but 
it gives the nonlinear Schrodinger equation. In practice, choosing a complex coefficient allows one to study 
the Hall Effect in superconductors, but we will generally choose F such that the time derivative of ip has no 
factor in front. The implication is that is only loosely associated with the density of superconducting 
electrons. What we have instead is a non-conserved order parameter. Our version of a continuity equation 
is 

which is more of a balance between energy in the system and dissipation in the current. 



1.4 Dimensionless Form of Landau-Ginzburg Equations 

It now falls to us to rewrite the Landau-Ginzburg equations in dimensionless form. In the process, we will 
expose the important length scales in the problem. The two length scales in the system are the magnetic 
field penetration depth and the superconducting coherence length. The illustrations of length scales from 
the time-dependent Landau-Ginzburg equations are very simple and are shown here much as they appear in 
Du, Gunzburger, and Peterson |]3^ . 

The magnetic penetration depth. A, is a measure of how far an applied magnetic field penetrates a sample. 
If we look at a perfectly conducting sample occupying a halfspace and apply a constant magnetic field, then 



the time-dependent Landau-Ginzburg equations, Eqns. 1.14 and 1.15, simplify to 



— VxVxA = VoA. (1.29) 



Taking the curl of this equation, we arrive at the London equation, 

1 

A2 

where the magnetic penetration depth is 



VxVxH+-'2H = (1.30) 



^ = \hr^- (1-31) 

The London equation predates the Landau-Ginzburg equations. It says that a magnetic field applied to the 
surface will have a decay length of A. 

We can make a similar calculation for the coherence length. If we imagine a metal normal to the left, 
superconducting to the right, then the order parameter has to rise from i/; = to V' = v— The 
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A = 



9 



1/2 



A / /3 mc 
~ ~ V 2TT~en 

H = V2HcH' 



a 



-9 



2 X 1/2 



2mQ! 



cH, 



' j' 



2V27rA 



47ra^ 



2\ 1/2 



a; = Ax' 
A = V2\HcA' 

t = T^t' 



Table 1.1: These are the rescalings of the time-dependent Landau-Ginzburg equations to their standard 
dimensionless form. 



Landau-Ginzburg equations in the absence of magnetic field reduce to 



2m 

The solution to this equation with the desired boundary conditions is 



— — tanh 



where 



-2ma 



V2^ 



1/2 



(1.32) 



(1.33) 



(1.34) 



We call ^ the coherence length. It is the order parameter's recov ery length to perturbations. 

Using the length scales shown above, we can (see Appendix A2) rescale the variables so that only two 
remain 



7 



dip 
'dt 



+ iK(f>ip + V - IV'I V - (- V + A)^xP = 



OA i 

— + V(I) + V xV X A- —{ipV^* - V'* VV-) + AIV-P = 0. 

The switch to dimensionless variables is nicely summarized in Du, Gunzburger, and Peterson 
complete summary is shown in Table |l.l| . The Onsager coefficient has become 



7 = 



aVtQ 



(1.35) 
(1.36) 
A more 

(1.37) 



where io is the rescaling of time shown in Table 1.1, and calculations from BCS theory determine it to be 
between 0.8 and 1.0 As you may guess, we'll take it to be unity. 
The other parameter is the Ginzburg-Landau parameter. 



A 



(1.38) 

which is a dimensionless parameter characterizing the metal under consideration. In this formulation of the 
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Landau-Ginzburg equations, it is the only such parameter. A metal with small k was once called hard — it 
resists penetration by magnetic fields.^ Metals with k < l/\/2 are called type-i, those with k > l/y/2 type-ii. 
In these new variables, the current is rescaled so that the nondimensional current 



7^(V'VV'* - V/VV^) - IV^^IA. (1.39) 

ZK 



The new boundary conditions are 



-Vip + A^Pj - n = (1.40) 

VxAxn = Hxn. (1-41) 

They still express the requirement that no supercurrent flow out of the sample and that the tangential 
magnetic field be continuous at the boundary. 



These nondimensional TDGL equations, Eqns. 1.35 and 1.36, are also invariant under a gauge transfor- 
mation, 

i/^^Ve'"^, A^A + Vx, and 0->,/,-^. (1.42) 

ot 

This indicates that, as written, the equations are not well-posed because they do not have a single set of 
solutions but a whole family of solutions joined by the function x- While it is possible to solve the equations 
analytically for a family of solutions without specifying a specific x, we will work with specific functions x, 
or specific gauges. We will derive those gauges as they apply in later chapters. 

We have introduced the Landau-Ginzburg equations from general properties of the free energy and with 
little reference to the physics of the electronic system of the normal or superconducting metals. The majority 
of this work will be concerned with the equations themselves. As just shown, the physical system informed 
the derivation at two points — the choice of SO{2) symmetry in the original free energy and the length scales 
of the gauge variables. A, ( and (j). We use those length scales to analyze the superheating of a half-space in 
the next chapter. 



''These days hard superconductors refer to type-ii superconductors with lots of pinning impurities to hold vortices. 



Chapter 2 

Superheating of Superconductors 



2.1 Introduction 

A first analysis of the free energy of a metal examines the properties of a bulk specimen large enough 
that surface tension associated with edges is negligible. Below a thermodynamic critical field, He, a bulk 
superconductor in equilibrium will be in a Meissner phase. Above that critical field, a type-i superconductor 
will become normal, allowing magnetic fields to penetrate, and a type-ii superconductor will enter the 
Abrikosov phase, characterized by the infusion of flux into the sample in the form of discrete vortices. 

If, instead of estimating the critical field for a phase transition just by comparing which phase has the 
least free energy, we, in addition, ask at what applied field flux at the boundary of the sample (with vacuum) 
can push past through the surface currents, we see that the phase transition is delayed. If one were to 
begin with a superconducting metal and steadily increase the magnetic fleld, there would not be a phase 
transition at He but for some higher value we call the superheating field, Hsh- (While we are not actually 
heating the sample, the word "superheating" refers to the similar delay of a phase transition in more familiar 
systems like water.) For applied magnetic fields well above the critical field, the bulk would prefer to change 
phases, but the order parameter at the surface decreases from its value in the bulk in order that the metal 
remain superconducting. On a microscopic scale, the depression of the order parameter is caused when the 
penetrating magnetic field aligns spins in the Cooper pairs. The phase transition is delayed only when the 
edges of the sample are considered. 

If we want to understand the time-dependent collapse of superconductivity in a superconductor, charac- 
terizing the superheated Meissner state will tell us our initial conditions for the collapse. The superheating 



solutions to our system (described in Eqns. 2.13-2.16) measure energy barriers separating local equilibrium 



from catastrophic collapse. Scaling solutions for cascading phase transitions often depend on first integrals 
of initial conditions during a superheated state. We will not be able to produce those scaling solutions, but 
a later chapter pursues a similar system found more tractable. 

The matched asymptotics in this chapter were published in a journal paper pTj . New material appears 
in the analysis of perturbations on the solutions. 



2.1.1 The Physical System 

The superheating field is more practically important for quite a few reasons. First, it is one of the most 
common ways to determine the Landau-Ginzburg parameter of a type-i metal [ p^ . Second, there have been 
proposals in the past ten years to use superconducting granules to detect weak electromagnetic elementary 
particles, including WIMPS and magnetic monopoles ||8^. Here, the granules act as bubble chambers, 
flipping from superconducting to normal when a particle strikes. A review is given in Barone 

Both measurements of the Landau-Ginzburg parameter and detectors study small spherical or ellipsoidal 
particles rather than large slabs of superconductor. Larger samples generally show almost no superheating 



11 
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Figure 2.1: A superconducting sphere expels applied magnetic fields. The increased field at the equator 
is called the demagnetizing field. To order X/R, the field at the equator is 3/2 the applied field. 



A [A] UA] 

Al 500 16000 0.03 
Sn 500 2300 0.2 
Pb 400 830 0.5 

Table 2.1: These are some sample coherence lengths and penetration depths for a few materials 
The dimensionless constant k = A/f . The coherence length and penetration depth appropriate to the 
Landau-Ginzburg equations differ slightly from the physical quantities. A good discussion of this is 
in Feder and McLachlan pal . 



due to defects. When early theories predicted almost no superheating Q, experiments on bulk supercon- 
ductors supported the results 0. It was then discovered that it is possible to make colloids of very 
uniform superconducting particles ranging in size from 10-100 /xm. Because the defects are smaller than the 
order parameter coherence length, superheating is near ideal. 

When calculating the superheating field of a spherical granule, one needs to account for the demag- 
netization fields caused by expulsion of flux from the granule. For a spherical grain which is perfectly 
superconducting, the field at the equator of a grain is three halves the applied magnetic field. Ellipsoidal 
grains and penetration depths complicate the exact calculation, but the experiment still rests on the physics 
of a field penetrating the surface of a superconductor. 

The master of the superheating field is Hugo Parr. He ascribes his lineage to Feder and McLachalan 
and in 1979, he published his final paper of a series on the subject (8|, ||, |6|, |8|||, |§. Using small grains 



of various metals, he found nearly ideal superheating fields in both type-i and type-ii superconducting grains 
in suspension. While he did not use SQUID techniques common to later studies on granules as detectors, 
his results showed that the prediction of the large superheating fields for type-i superconductors is accurate 
to within a percent. He used the granule system to study the temperature dependence of the penetration 
depth and coherence length in detail. He even noted the recently rediscovered "giant vortex state" of type-ii 
granules not large enough to contain a single discrete vortex Q . 

We want to calculate the superheating field, i?sh, as a function of k, the Landau-Ginzburg parameter. 
That parameter is the ratio of the coherence length of the Cooper pairs, ^, to the magnetic field penetration 



depth, A. This is known to within 10 % for most superconductors. A few are shown in Table 2.1. We will 
examine the most basic system we can, a superconducting half space. 
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per unit 
length 




Nucleus Radius 




Figure 2.2: A small nucleus of the normal state in a superconductor will have a large positive surface 
tension. No nucleus can grow unless its size is larger than a critical radius where it is energetically 
favorable. 



2.1.2 A Rough Estimate 

If we focus on the first entry of flux into a superconductor, we can think of a smaU nucleus at the edge of the 
metal. For a type-ii superconductor, the flux will enter as vortices typically much smaller than the penetration 
depth of the magnetic field. In 1964, Bean and Livingston estimated the critical magnetic field at which 
flux would enter a type-ii superconductor by imagining a small test flux interacting with the unperturbed 
system The same picture of a small flux line interacting with a nearly intact superconducting state is 
still used to refine magnetization curves |jl^. For a type-i superconductor, we can model the entry of flux 
as the nucleus of a phase transition from superconducting to normal. When the critical radius of a nucleus 
of normal metal becomes as small as the typical perturbation of the system, a phase transition will occur. 
Our goal is to determine the transition field as a function of the Landau-Ginzburg parameter, k. 

Picture a nucleus of normal metal in a superconducting cylinder. We can look at a two-dimensional 
slice of an infinite cylinder so that the nucleus has an energy per unit length determined by the competition 
between volume energy and surface tension. Surface energy per unit volume depends on the difference 
between the applied magnetic field and the critical field 



The applied field is assumed to be larger than the critical field for superheating. The energy cost of an 
interface between normal and superconducting regions is the surface tension. We can write a nondimensional 

form of the surface tension 



unit volume 



E 




(2.1) 




(2.2) 



where Osborn and Dorsey [§2| determined that, for asymptotically small k, ans is 



0.943 



0.880 



(2.3) 



K 



From the competition between the surface and bulk energies, we can write the total energy of a nucleus of 
normal metal 
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The radius at which this energy becomes negative is the critical radius for nucleus growth. 

A graph of the energy dependence is shown in Fig. 2.2. Smaller nuclei have a large ratio of surface to 
volume, so the positive surface tension makes them shrink. Nuclei large enough for the negative bulk energy 
to become dominant have radii larger than the critical radius 

We can estimate that a phase transition must occur when the critical radius of the nucleus is about the same 
size as the penetration depth 

- A (2.6) 

Solving this for the applied field where the system must collapse, we find 

= H^{l + 2ans). (2.7) 
If we insert our asymptotic value for ans, we find 

1 34 

Ho/H,^ —j^~OMl. (2.8) 

The factor in front is about twice the value calculated later, but the order of the first term is correct. The 
superheating field of a type-i superconductor will diverge as k — s- with a power of k~'^/^ because the positive 
surface tension diverges as k~^. 

Now we look at a more mathematical approach which accounts for surface currents and magnetic fields 
using the Landau-Ginzburg equations. 



2.1.3 Formulation of the Problem 

The bare requirements for a superheating problem are an applied magnetic field and a surface. We will study 
a superconducting half-space where all space a; > is filled with metal and all space a; < is empty with an 
applied magnetic field H. This is not a time-dependent problem so we need use only the time-independent 
Landau-Ginzburg free energy 



G 



1 



d-'x i -iVr + 7:\W + (V X A)' - 2V X A • H - 



1, 



- V + iA V 



(2.9) 



where k is the GL parameter, tjj is the amplitude of the superconducting order parameter, A is the vector 
potential (B = V x A), and H is the magnetic field applied to the side of the sample. The lengths are in 
units of the penetration depth A and fields are in units of \/2Hc- While supercurrents shield the metal from 
applied magnetic fields, they are not dissipative so we don't need to include any terms related to the electric 
field in order to model this system. 

The nondimensional free energy is invariant under a local linear transformation of the form 



and 



(2.10) 



The gauge transformation amounts to an extra degree of freedom. It is called local because the transformation 
depends on coordinates, x = x(^)- While we could, in principle, perform calculations without specifying 
a gauge, the easy and common practice is to fix the gauge. We fix the gauge by completely specifying a 
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function x in terms of observables. In this case, if one writes tp = /e'^, we can make "0 real by defining 

X^^e/K ^ S,=iJe^' ^fe^'e-^' = f. (2.11) 
The phase of tp is ehminated without any loss of generality. The transformed free energy is 

G = y d^x I -1 + Q'/' - + i./^ + ( V X Q)2 - 2( V X Q) . h| . (2.12) 

where / = ji/'l is the amplitude of the superconducting order parameter and Q is the gauge-invariant 
vector potential. Having transformed the free energy, we take the functional derivative in order to get the 
equilibrium equation 

\v^f-fQ^ + f-f = Q (2.13) 



VxVxQ + /"Q = (2.14) 



with boundary conditions 



V/-n = (2.15) 
(VxQ-H)xn = 0. (2.16) 

If we think of a classical mechanical system where x —> these would be equations of motion for a particle 



in a potential. Eq. 2.14 describes the currents in the system in the form 

Jtotal Jsupcrconducting — (2T7) 



because there is no normal current. Furthermore, the divergence of Eq. 2.14 demonstrates explicitly that 
the supercurrent is zero. 

Now we further restrict our treatment to a single dimension by treating independent variables as constant 
in directions parallel to the surface. Since our sample is a halfspace a; > 0, we apply a magnetic field parallel to 
the surface by specifying H = haZ. We represent that magnetic field with a vector potential Q = (0, q{x),Q). 
All spatial variations are in the x direction so that the resulting equations are 

\f" - q\f + .f- f = 0, (2.18) 

q" - fq - 0, (2.19) 
h = q'. (2.20) 

where primes denote derivatives with respect to x. These much simplified equations will suffice us the whole 
of the derivation of the superheating until we need to examine its stability. 

2.1.4 Previous Work 

The Landau-Ginzburg equations have been used to study the superheating field since their inception. All 
work done with them has been approximate because these are coupled nonlinear equations and intractable 
for this application. All approaches focused either on large or small asymptotics in n. We will look at both, 
discussing in greater detail the work on small n. 

In 1950, Ginzburg and Landau wrote Eqns. 2.18 and 2. IS in order to find the effect of a magnetic 



field applied parallel to a superconductor. Assuming the order parameter varied only a little for a hard 
superconductor (small k), they wrote / = / + ./ where f — I was the bulk solution and / < showed 
the depression of superconductivity at the boundary. Second-order terms, and qf, were negligible. The 
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resulting equations were 



^/"- 2/ -9^ = (2.21) 
<?" = q- (2.22) 
Solving these is straightforward, and the approximate order parameter is 

While they did not further investigate the superheating field, the above equation did establish that the order 
parameter at the surface responds to applied magnetic fields according to kH^. Note also that our dimen- 
sionless units measure distance in terms of the magnetic penetration depth so that this first approximation 
shows q oc e~^. 

In 1958, Ginzburg [Q estimated the superheating field for large and small k. For k — > oo, he solved the 
Landau-Ginzburg equations exactly to find Hgh ~ l/v^ (Recall the field is normalized to V^Hc). He did 
not examine the stability of that solution to discover the physical superheating field is actually lower. 

For small n, however, Ginzburg made fruitful conjectures from scaling arguments: "From [Eq. ^.23| ], and 
from an analysis of [Eqs. ^.18 and 2.19| with the introduction of variables C — ^/kx, x = f /V^ ^^'^ ^ = v^'i'i 



it can be inferred that, for ^/K <C 1, Hsh — const/y^." How do we explain these choices? 

We want to isolate and solve the relevant parts of the coupled nonlinear differential equations. Here, the 
relevant parts are those concerning the magnetic field. For hard superconductors, the order parameter will 
remain nearly constant while the magnetic field dies off quickly. We also have the sense that the superheating 
field for a hard superconductor is large. A good first step would be to rescale the vector potential so that it 



is closer to being 0(1). Given choices like b ^ Kq or b ~ \/l^q, one might look at kH^ in Eqn. 2.23 and choose 
the latter so that kH^ would become the 0(1) (V x b)^ in the new variables. There is also a mathematical 
expression for our expectation that the order parameter remain steady while the vector potential varies. We 
simply change the length scale of the problem so that the vector potential looks more like a perturbation. 
With units in terms of the penetration depth, we have q oc . We can use units C = \/kx so that C remains 
0(1) when x is large. We now have the following set of equations 

-r --b^f + f-f = (2.24) 

K K 

Kb" - fb ^ 0. (2.25) 

Since we claim it is the behavior of the vector potential that describes this problem, we want to include the 
second equation. We can do that by rescaling y^X = / with the effect 

-X" --b^ + X-^x"" = (2.26) 

K K 

b" ~ x^b = 0. (2.27) 



From here, one can do series solutions for x b in k. The zeroth order solution is that x = 1/v^ 



Just from the zeroth order solution, we can look at Eqn. 2.27 to find b — —Ce^^^'^ — — Ce^^ so that, 
q = bj ^/k = —{C I ^fK)e~'^ . Our rescalings along the way sought to ensure that all proportionality constants 
were 0(1). We conclude, therefore, that H — const/y^. Ginzburg determined that constant numerically to 
find 

H-sh — — i^Hc- (2.28) 



The Orsay Group on Superconductivity |Q (including de Gennes), in 1966, took a similar physical 
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approach with a different mathematical technique. Rather than look at scaling, they used a variational 
argument. The basis is that the magnetic field goes rapidly to zero while the order parameter takes its time. 
The order parameter is calculated in the absence of the applied field 

4/" = /(/' - !)■ (2-29) 

The vector potential is calculated for a fixed order parameter so that 

q ^ g(0)e--^'(")". (2.30) 

Then they put these estimates into the free energy and concentrate, as did Ginzburg, on the terms involving 
q, q' , and /'. Terms in / are dropped. They vary the free energy with respect to /(O) to find what /(O) 
gives a free energy minimum. The result is a nontrivial estimation of the applied field in terms of /(O) 

= — /^0)(l-/2(0)). (2.31) 

K 

The applied field always reaches its maximum for /(O) = 1/V2 where 

H,h = 2-1/4^-1/2^,. (2.32) 
This exact analytical coefhcient is the same as that derived later in this chapter, but the eponymous Ginzburg 



would cite our paper from 199612^ in his 1998 work on the thermoelectric effect We can see in Eqn. 2.31 
(plotted in Fig. 2.11 on page ^5[) that the superheating field occurs where — 0. It is also apparent that 



the order parameter at the surface is nonzero at the maximum applied field. 

More recently, Hugo Parr combined intuitive rescalings and a variational approach to derive the next 
term in the superheating field. 

i/., = 2-V4^-V2 + (2.33) 

The techniques are not significantly different from those above. The reasons for his rescaling are not clear 
and the calculation for the second term is very complicated. His work is important because the addition of 
the second term brings the asymptotics quite close to the full solution of the Landau- Ginzburg equations 
as shown in Fig. on page |2^. This permitted him to make practical verifications of the estimated 
superheating field. 

Our contribution to the examination of the superheating field is two-fold. First, we determine the 
superheating field to greater precision than previous work. Second, and more importantly, we use a method 
called matched asymptotics to do that derivation with surety. While matched asymptotics are informed by 
the physical considerations discussed above — short penetration depth, slow variation of order parameter — the 
method transforms the problem into linear series solutions which can be integrated mechanically. 



2.1.5 Boundary Layer Method for a Converging Channel 

The boundary layer method is a standard technique of singular perturbation theory, but it is much less well 
known than WKB, so we will review it here. An excellent book on the subject is Bender and Orszag [ p^ . 
For physicists, a more amusing introduction might be a simpler problem from fluid dynamics. 

Landau and Lifshitz fz^ derive the exact solution of the equations of motion for a viscous, incompressible 
solution flowing into a converging channel, but they don't find a closed form solution. We follow their brief 
derivation enough to make sense of the main differential equation, then use the boundary layer method to 
find a closed form solution for low viscosity, high Reynolds number fluids. This is a system with a single 
small length scale, identified by the Reynolds number. It should be good practice for the Landau-Ginzburg 
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equations which exhibit two or more length scales in multiple coupled nonlinear equations. 

The system geometry is shown in Fig. 2.3. The fluid flows into the channel and out a pinhole. We assume 
the flow is entirely radial. We could determine the equation of motion for an inviscid liquid just from mass 
conservation. Given a mass Q of liquid per unit time, 



Q = P 



a/2 



a/2 



vr d(j). 



For an inviscid liquid, the solution would be 



-pvr 



2ttQ 



2tt par 
In order to include viscosity, we will use radial coordinates where the Navier-Stokes equations are 



dv 1 dp 

= ^ + " 

or p or 



d^v 1 d^v 1 dv V 

Qj,2 j,2 g^2 J, Qj. j,2 

1 dp 2v dv 



pr d(p r 
d{rv 







dr 



= 0. 



(2.34) 

(2.35) 

(2.36) 
(2.37) 
(2.38) 



These equatio ns a re still supplemented by the mass conservation equation, Eq. |2.34| . Because of the last 
equation, Eq. 2.38 , we know v cx l/r, so we change variables to 



rv 



Substituting u in Eq. 2.37, we can integrate to flnd the pressure 

121/2 



1 
P 



[P - Pwall) 



(2.39) 



(2.40) 



Substituting the pressure into Eq. 2.3£, we flnd 
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Figure 2.4: Fluid flow is uniform in the center of the channel. It drops to zero only near the sides. 
The region affected by viscosity is C'(l/\/Re). The estimated stream velocity for a large Reynolds 
number, Re = 43, is a very good approximation of the exact solution. 



Noticing that the left-hand side depends only on (p and the right-hand side only on r, we assign each to a 
constant, 2c, and arrive at the main differential equation for this system 

u" + 4m + 6u^ = 2c. (2.42) 

Here, u depends only on (f) and c is a constant of integration. 



Landau and Lifshitz express Eq. |2.42| in integral form and discuss the nature of the solutions. We will 
examine the case where the viscosity, v, is small, or Reynolds number, Re = IQl/pv, is large. Here we 
expect that, through most of the channel, the liquid will move with nearly the same velocity as the inviscid 



case, Eq. 2.35 . Near the sides, however, the velocity at the walls must drop to zero. The affected region in 
the channel will be ©(l/v^Re). The affected region is our boundary layer. We can rewrite the differential 
equation in different variables to see why there is a boundary layer. 

First, note that the constant, c is 0(Re^). This could be proven by estimating Pwaii from Bernoulli's 



Equation. We also know from Eq. 2.35 that the maximum u is less than 



=-Re-^. (2.43) 



ivpar 6ar 

That tells us that u is C'(Re). Let's first write the main differential equation in terms of 0(1) variables. 
Substitute U = -u/Re and C = c/Re^ to find 

^U" + —U ^W"^ = 2C. (2.44) 
Re Re 

It is customary to write the small parameter as e = 1/Re 

dJ" + Aeu + = 2C. (2.45) 

While this equation applies in the center of the channel, it is our outer equation because it describes the 
region far from the boundary layer. For flows near the center, U" is small so we don't expect it to be 
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important. We will need to rescale in order to look at flows near the side. We can use a simple series 
solution to solve the outer equation. Substituting the series 

U = UQ + eUi + e^U2 + ■■■ (2.46) 

C = Co + eCi+e^C2 + .... (2.47) 

we find 

6U^ = 2Co. (2.48) 



The solution is Uq = \/ Co/3 in the center. We can write down the next order differential equation, too 

f/^' + 4{/o + 12;7o{/i = 2Ci ^ f/i^^-i. (2.49) 

0(70 'J 

We determine constants, Co and C\ from matching the solution near the boundary. A nice side effect of the 
series solutions is that each succeeding order is always a linear differential equation. 

Near the boundary, we expect the second derivative of V to be important. We express this by rescaling 
the (/)-coordinate, <I> = 4>/^. 



AeU + W^=C. (2.50) 



We have the option of rescaling U , in addition. While we expect U is smaller near the boundary, rescaling 
U would give us U" = C which could not fulfill boundary conditions. We usually seek rescalings 

where more than one term is of the highest order of e because you usually need two or more terms to find a 
non-trivial solution. Scaling variables such that nontrivial solutions exist is called the principle of dominant 
balance. Currently U" is balanced against QU^ — C which seems like it could fulfill boundary conditions 
non-trivially. If we substitute series into this equation, we find 



U = C/„o(-2 + 3tanh2(±v/-3C/o($ - k)). (2.51) 

The integration constant, fc, is determined from matching to the outer solution. 

Now that we have an inner and outer solution, we match them. We do that by writing the inner solution 
in terms of the outer variables, <I> = (t>/\/e. 



U = C/,„o(-2 + 3tanh^(±v/-3C/o(e-i/2(/) - k)). (2.52) 

Now we define the constant k such that [/(—a/2) = 0. The inner equation is already set to match J/outor = 
UmQ. We therefore see the real size of our boundary layer is l/y%5|[/o| — 1/ •\/3|uo|R.e. 

The solution for the boundary layer meets the boundary condition, u = 0, at the side of the channel and 
meets the solution to the inner equation as — *■ oo. The outer solution would not, by itself, be capable of 
fulfilling the boundary conditions at the walls, but it does match to the inner solution near the walls. The 
solutions of the superconducting half-space will have the same property. Inner and outer solutions together 
form a uniform solution to the system. 



2.2 Asymptotic Expansions for Small-AC 

In this section we will develop an asymptotic expa nsion for the superheating field for small-K, using the 
method of matched asymptotic expansions jl^, |l05|l . For small-K the longer length scale is the coherence 
length f . In order to treat the problem with singular perturbation theory, we need to write it such that 
the factor in front of the highest order derivative is small. Second order equations with a small parameter 
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Figure 2.5: The estimated stream velocity for moderately large Reynolds number, Re = 9, is a fair 
approximation of the exact solution. 



modifying the highest derivative tend to exhibit rapid behavior near the boundary. In this case, the rapid 
behavior is the penetration of the magnetic field. That field is the boundary of boundary layer theory, which 
we use here. 

Rescaling x by k, we introduce a new dimensionless coordinate x' = kx which means we now measure 
distance using the coherence length. The resulting GL equations in these "outer variables" are 



/" - + ,/-/' = 0, 

h — K,q' , 

with the primes now denoting differentiations with respect to x' . 

Outer solution. In order to obtain the outer solutions expand /, g, and h in powers of k: 

f = k + ^^h + >i\f2 + . . . , 

q = qo + nqi + K^q2 + . . . , 
h ^ ho + nhi + K^h2 + • . . - 



Substituting into Eqs. ( 2.53 )-( 2.55 ), at 0(1) we have 

/o - qlfo + fo- ,/o' = 0, 
-foqo = 0. 



(2.53) 
(2.54) 
(2.55) 



(2.56) 
(2.57) 
(2.58) 



(2.59) 
(2.60) 



Since we want / 1 a,s x' ^ oo, the only possible solution to Eq. ( 2.60 ) is go — 0. We can then immediately 
integrate Eq. ( 2.59 ), 



fo{x') = tanh 



Xo 



V2 



(2.61) 
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with xq — xo{n) an arbitrary constant.]^ To 0(k), the outer equations are 



/{' - 2qofoqi - qlh + fi ~ S/^/i = 0, 
-/o\i - 2/ogo/i = 0, 
ho = 0. 



(2.62) 
(2.63) 
(2.64) 



Once again, the only solution to Eq. (2.63) is qi — 0; substituting this into Eq. ( 2.62] ), we find /i = Ci/g, 
with C'l a constant: 



(2.65) 



(n) 

We can continue in this manner; at every order q^ — 0, hn = 0, and /„ ~ Cnf^ , with the C„'s constants 
which are determined by matching onto the inner solution. 

Inner solution. The outer solution breaks down within a boundary layer of 0(k) near the surface. This 
suggests introducing a rescaled inner coordinate X — x' / n, so that X = 0{\) within the boundary layer. 
It is also possible to rescale / and q, with the hope that this will lead to a tractable inner problem. Such 
a rescaling must lead to a successful matching of the inner and outer solutions; i.e., the inner solutions 
as X — > oo must match onto the outer solutions as x' — s- 0. Since /o(0) — tanh(2;o/\/2), then assuming 
that a;o 7^ we have /o(0) = 0(1), indicating that the order parameter should not be rescaled in the 
inner region; therefore we set f{x') = F{X) in the inner region. However, from the outer solution for the 
vector potential we see that the only constraint on q{X) in the inner region is that q{X) — )■ as X — )■ cxi 
(presumably exponentially). Therefore, we are free to rescale g by k in the inner region, hopefully in a way 
which simplifies the inner equations. One possibility is q{x') — k^°'Q{X)] substituting this into the GL 
equations, Eqs. (2.53)-(2.55), we see that unless 2a is an integer, fractional powers of k will be introduced 
into the inner equations, contradicting our expansion of / and q in integer powers of k in the outer region. 
Therefore, the simplest assumption is that a — 1/2, leading to the following choice for the inner variables: 



x' = kX, f{x')^F{X), q{x')^ K~^/'^Q{X), h{x') = H{X). 

In these variables Eqs. ( 2.53D -( 2.55| ) become 

F" - kQ'^F + k2(F - F^) = 0, 
Q" - F^Q = 0, 
n^'^H = g', 

where now the primes denote differentiation with respect to X. The boundary conditions are 

F'{0) = 0, H{{)) = Ha. 
The next step is to expand the inner solutions in powers of k: 

F = Fa + KFi + K^F2 + . . . , 

Q = Qo + i^Qi + K^Q2 + • ■ • , 

H = K^^/^Ho + K^'^Hi + ... . 



(2.66) 

(2.67) 
(2.68) 
(2.69) 

(2.70) 

(2.71) 
(2.72) 
(2.73) 



Note that there is no term of 0(1) in the expansion for iJ, since we would be unable to match such a term 

^Because xq depends on k, it should also be expanded inaeries. For the orders of terms used in the following derivation, the 
differences are negligible. A letter from the BoUev et al. |13| pointed out the difficulty which was painstakingly recalculated 



and corr ecte d by Joh 



in Table 2.2 and Eq. 2.94 



X2i Bartolo in an erratum |p8(| . The affected terms are of high enough order that they show up here only 
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to the outer solution. Using the boundary condition H{0) = Ha leads to 

Ha = K-^/^HoiO) + K^/^Hi{0) + ... . 
Substituting these expansions into Eqs. (2.67)-(|]69), at 0(1) we obtain 

Fq' = 0, Qo — F^Qo = 0, Hq = Qq. 



Solving these equations subject to the boundary conditions (2.70) (we also need Qq 
to match onto the outer solution), we obtain 



as X 



Fo{X) = Ao, Qo{X) = Boe 



-AnX 



Ho(0) = -AoBo, 



(2.74) 

(2.75) 
oo in order 

(2.76) 



with Aq and Bq constants. In what follows we will assign i^„(0) — An and (3„(0) = B„ for notational 
simplicity. The 0{k) equations are 



Pi — QqPo, Q'i — FqQi = 2FoQoFi, Hi = Q\. 



Solving with the boundary condition i^i'(O) = 0, we obtain 

Fi{X)^Ai 



-2AoX 



1]. 



-AoX 



Bi 



16^2 



-2AoX 



-16 



i7i(0) 



AlAi 



IBl 

8 An 



2 X + AAlX^ 



AqBi - AiBo. 



Finally, to O(k^) we have for F2 



-Fo + F^ + 2QoQiFo + QlFi, 



(2.77) 
(2.78) 

(2.79) 
(2.80) 

(2.81) 



the solution of which (with i^2(0) = 0) is 
F2{X) = 



YI_Bl_ 
128 Al 



1 BlAi 
4 Al 



2 Ao 



A2 



IBoBi IBlAi 
2 ^0 



4 Al 



AM 
32 Al 



IBt 
8AI 



1 BlAi 

2 ^0 



X--Ap{l-Al)X^ 
8 Ao 



-4AoX 

128 Al 



(2.82) 



The expression for Q2 is even more unwieldy, and is not needed in what follows. 

Matching. To determine the various integration constants which have been introduced we must match 
the inner solution to the outer solution. Since the outer solution for q is simply q = 0, and all of our inner 
solutions decay exponentially for large X, the matching is automatically satisfied for q, as well as for the 
magnetic field h. To match the inner and outer solutions for the order parameter, we are guided by the van 
Dyke matching principle [105|, which states that the m term inner expansion of the n term outer solution 
should match onto the n term outer expansion of the m term inner solution. While the matching principle 
works for any pair, (m, n), experience shows which choices yield meaningful results at a particular order of 
the computation. In our case we will take m = 3 and n — 2. Therefore, write the two term outer solution 
/o(x') + Kfi[x') in terms of the inner variable X, and expand for small k, keeping the first three terms in 
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the expansion in powers of k: 

foinX) + nfiinX) 



tanh 



V2 



+ K scch 



2 ^0 



V2J V2 



-K^sech^ I 

.V2 



tanh 



V2 



[Ci + X] 



-CiX- 



X' 



(2.83) 



Next, write the three term inner solution Fq{X) + kFi{X) + k^F2{X) in terms of the outer variable x' , and 
expand for small k, this time keeping the first two terms of the expansion: 



Fo{x'/k) + kFi{x'/k) + k^F2{x'/k) 



Ao + ^x'-IAo{1^AI)x'' 



(2.84) 



By writing both expressions in terms of x' , and equating the various coefficients of x' and k, we see that the 
expansions do indeed match if we choose 



^0 = tanh 



V2 



Ai = -f- + sech2 
4Ao 



Bo = -2I/W (^-^^ = ~2i/4(l ~ AD' 
^/2y ^/2 4 



/2 



i2n C*! 



B - ^ ^" V2A„{1-Al) Ci 



Eliminating C'l 



V2A0A1 , 3 53 1 1 _ ^; 

-Dl = H — -To + - 



Bq 32 Aq 2 i?o 
Substituting into our expressions for Hq{0) and i?i(0) from Eqs. ( ^.7(: ) and ( ^.80 ), we obtain 

i^o(0) = 2l/4^o(l-^g)^/^ 

2^/4 (2A2 + 14)(1-A2)V2 2^/4(2^2 

-Wl((j) = ^7-; ^ 7; ^1- 



64 



^0 



(1-A2)V2 



(2.85) 
(2.86) 
(2.87) 
(2.88) 

(2.89) 

(2.90) 
(2.91) 



In order to calculate the superheating field (or, more correctly, the maximum superheating field), we need 
to maximize -ffo(O) and i?i(0) with respect to Aq and Ai. Maximizing Hq{0) with respect to Aq, we find 
that the maximum occurs at Aq = 1/V2, Bq ~ — 2~^/^, so that i?o(0) = 2~'^/''. Substituting this result 
into -ffi(O), we find the surprising result that the coefficient of Ai is zero, and Hi{0) = 2^/415/64. Our 
superheating field is then 



2-3/4^-1/2 



(2.92) 



In order to determine Ai we need to proceed to a higher order calculation. The method is the same as 
before, although the algebra quickly becomes tedious; we have used the computer algebra systems Maple V 
and Mathematica to organize the expansion. The results from a six term inner expansion are summarized 
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_n An Bn Xn H,, (0) 

"O 2=^ ^F^74 2i/2tanh-^(2-i/2) 2^3/3 

1 -7/32 -(9/16)21/4 -(15/16)21/2 (15/64)2^/4 

2 (395/2048)21/2 (147/512) 429/512 -(325/2048)21/" 

3 ? ? ? (14191/65539)23/4 

4 ? ? ? -(67453267/62914560)21/4 

Table 2.2: The lower order integration constants are shown in the text. Higher orders required 
symbolic mathematics software. Miraculous cancellations allowed derivation of the first six terms of 
the superheating field without knowledge of most third- and fourth-order constants. 



in Table 2.2. Including the next order term in the expansion in the superheating field, we have 

15V2 325 



H^^ = 2-3/4«-i/2 



1 



32 



1024 



k2 + 0(k-^) 



(2.93) 



The first term is exactly the result obtained by the Orsay group ^ , who used a variational argument 
to obtain their result. The second term is identical to the result obtained by Parr Parr combined 

an inspired guess for the behavior of the order parameter near the surface with a variational calculation in 
order to obtain his result. It is interesting to note that our result for Ai also agrees with Parr's result. The 
advantage of the method of matched asymptotic expansions is that we can make this expansion syst ematic, 
and therefore in principle carry out this expansion as far as we wish. The third term in Eq. ( |2.93D IS one 
of the new results of this chapter; the fourth and fifth terms are included in Table 2.2. With the five-term 
expansion for iJgh it is possible to employ resummation techniques to improve the expansion. For instance, 
the [2, 2] Pade approximant |l^ is 



1 + 4.6825120 K + 3.3478315 
1 + 4.0195994 K + 1.0005712 k2 ' 



(2.94) 



In Fig. 2.6 we compare the numerically calculated superheating field against the one, two, and three term 
asymptotic expansions. The one term (i.e., the Orsay group) result never seems particularly accurate. There 
is a marked improvement with the two term expansion, with the three term expansion offering only a modest 
additional improvement. The [2, 2] Pade approximant agrees with the numerical data to within about 1 % 
all the way to k = 1. 

Uniform solutions. From the inner and outer expansions it is possible to construct uniform solutions, 
which are asymptotically correct for all cc as k — *■ 0. To do this we simply add the inner and outer solutions 
of a given order, which guarantees the correct behavior in the outer region as well as in the boundary layer. 
However, this would produce a result which was 2/niatch in the matching region, so we need to subtract 
/match ill order to obtain the correct behavior in this region. As an example, we will construct the 2- 
term uniform solution for the order parameter. Adding the two-term outer solution, fo{x') + Kfi{x'), to 
the two-term inner solution, Fo{X) + kFi{X), subtracting the solution in the matching region, which is 
I/V2 + {^/2/4)kX — (15/32)fi;, and writing the entire combination in terms of the original variable x (which 
is the same as X), we obtain 



/unif,2(a;) = tanh 



xo 



V2 



15 
— f 
16 



-K sech 



- Xo 



V2 



4 



(2.95) 



As 2: 00, /unit, 2(2;) 1; also, /unif,2(0) = l/\/2 - (7/32)k, as we expect. However, /u„if^2(0) = (15/64)^2, 
so that the zero-derivative boundary condition is only satisfied to 0{k). 

In Fig. 2.7 and Fig. 2.S we compare the numerically calculated order parameter and magnetic field with 
the two term outer solutions and the three term inner solutions. The agreement is quite good for k = 0.1, 
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Figure 2.6: A comparison of the numerically calculated superheating field ffsh (heavy line) with the 
three term asymptotic expansion for small-zt, and the [2, 2] Fade approximant. The one-term expansion 
due to the Orsay group deviates systematically from the calculated superheating field. The two- and 
three-term expansions provide a marked improvement over the one-term expansion. 
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0.0 1.0 2.0 3.0 0.0 2.0 4.0 6.0 

X, outer length scale x, inner length scale 

Figure 2.7: A comparison of the three term inner and outer solutions for the order parameter and 
the magnetic field with the numerical solution for k = 0.1. The asymptotic solutions approximate 
the computed values only in the appropriate regions. The matching region where the inner and outer 
meet is 0{k) as can be estimated from the inner solution for /. 
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1.0 



m 

E 
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to 

Q. 
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0.0 



numerical 






"■-^ inner 
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2.0 



1.6 - 



0.0 1.0 2.0 3.0 
X, outer length scale 




0.0 2.0 4.0 6.0 
X, inner length scale 



Figure 2.8: The same as Fig. 2.7 for ft = 0.5 



with deviations appearing at k = 0.5. These figures also illustrate the existence of a matching region where 
the inner and outer solutions overlap; this region grows as k — > 0. Lastly, we show in Fig. 2.9 how the two 
term uniform expansion constructed earlier supplies a uniform approximation to the order parameter and 
magnetic field over the whole region. 



2.3 Stability Analysis of the Solutions 
2.3.1 Derivation of Equations 

We test the stability of the solutions by solving the second variation of the free energy to see whether 
the Landau- Ginzburg solutions sit in a minimum, maximum or saddle point of the free energy. The second 
variation of the Free Energy can be found any one of a number of ways. The first solution shown is equivalent 
to using Calculus of Variations and the second is known as a time-dependent formulation. 

If we perturb the extremal solution (/, q) of the GL equations by allowing /(x) — > /(x) -I- /(a;) and 
q{x) q(x) + q{x), then the second variation of the free energy functional is 



dx 



\P + (Sf + q'- 1)P + 4fqfq + + q 



~2 



(2.96) 



The boundary conditions on / and q should be chosen so as to not perturb / and h at the surface, so that 



/'(o) = q'{0) = 0, /V) - = 0. 



(2.97) 
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Figure 2.9: A comparison of the two-term uniform solution for the order parameter, /unif 2{^) (dashed 
line), witli tlic numerical solution (solid line) at k = 0.5. The disagreement of the uniform solution 
with the boundary condition at a: = is of order k'^. 
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We can then integrate Eq. (2.96) by parts to obtain 

1 



S^T^ I dx 





(2.98) 



This quadratic form can be conveniently written as 







5^T^ I dx{}.,q)Li ii\ (2.99) 



where Li is the self-adjoint linear operator 



n= ' 2-'^, K (2.100) 



9 / V 2/g 
In order to analyze the stability, expand / and q as 




, /^C„ i" , (2.101) 

where the c„'s are real constants, and (/„, (j„) is a normalized eigenfunction of Li with eigenvalue i?„: 

Li( ^.")^E,J ^J'). (2.102) 



9n / V * 



Then 



S^J^ = J2Encl. (2.103) 



The second variation S^J-' ceases to be positive-definite when the lowest eigenvalue first becomes negative, 
indicating that the corresponding solutions (/, q) of the GL equations are unstable. Therefore the entire issue 
of the stability of the solutions has been reduced to finding the eigenvalue spectrum of the linear stability 
operator Li, which in the k limit can be studied using matched asymptotic expansions. 

The condition that the second variation be positive definite is identical to the traditional requirements 
of stability calculations using Calculus of Variations. It measures stability of the system with respect to 
infinitesimal perturbations. We also know, cither from theorems of Calculus of Variations or from equivalence 
with the Schrodinger equation, that the solutions with lowest eigenvalues will never cross zero — they will be 
lowest energy bound states. 

Another way to analyze the stability of solutions to the GL equations is to substitute perturbations into 
the time-dependent GL equations and find whether they remain bounded over time. This is often presented 
as a second way to measure the same stability as that of the Calculus of Variations, but it is essentially 
different. Time-dependent stability analysis is dynamical. Let's look at the derivation of equations to see 
the difference. 

In many systems, one can determine the stability of the solutions by substituting time-dependent per- 
turbations 

fix) ^ /(x) + /(x)e-^* (2.104) 
q{x) q{x)+q{x)e-^* (2.105) 



into the time-dependent equations. Here (/, q) are the exact solutions. Then linearize the resulting equations 
in order to find linear stability. As above, positive eigenvalue, E, indicates a stable system. 
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In the gauge where ip is real, the set of time-dependent differential equations which would yield the exact 
same results as the Calculus of Variations approach are 



dl 

dt 



^ + Vx VxQ+/2q = 0. 
ot 



(2.106) 
(2.107) 



We just inserted time-dependence of the correct sign to look like a diffusion equation. If we put the pertur- 
bations into the above equations (Eqns. p. 104 and 2.105), the resulting equations are 



1 



r + {'Sf + q^-l)f + 2fqq = Ef 
-q" + fq + 2fqf - Eq. 



(2.108) 
(2.109) 



You see that the time dependence cancels when we linearize in / and q leaving the eigenvalue E. When 
i? > 0, the solutions (/, q) are stable with respect to infinitesimal perturbations. These equations are formally 



identical to Eqns. 2.102 



The Eqs. 2.107 are almost, but not quite, the full time-dependent Landau- Ginzburg equations. We 
haven't yet seen the TDGL in the gauge where "0 is real, so let's write them now. If we make a gauge 
transformation to Eqs. 1.35 and 1.3(; of the form 



Q = A + Vx 



dx 
dt 



where we specify 



KX 



for rp = fe 



i9 



we find the TDGL in the gauge where tp is real 







9Q „ 



7«:20/ + V • (Q^T/i) = 
+ V X V X Q + /2q - V X H = 0. 

Given our initial solution in one dimension is of the form 

Q = iO,qy{x),0) 
f = fix) 



we know the divergence on the right of the second equation, Eq. 2.113 

2/V • Q = Q • V/ - 



IS zero, 



(2.110) 
(2.111) 

(2.112) 
(2.113) 
(2.114) 



(2.115) 
(2.116) 



(2.117) 



which makes 4> = 0. The full TDGL then reduce to Eqs. |2.107| . This is true of our initial solution, but not 
true for general one-dimensional syste ms. 



In two dimensions, however, Eqns. p. 112 
degree of freedom in 



-2.114 



cannot be reduced to Eqn. 2.107 , They include an extra 
This extra degree of freedom allows them to include normal currents, which are 
dissipative. Recall that the Calculus of Variations perturbations do not include dissipation. Because the 
equations in one dimension are the same as the TDGL, they do describe dynamic perturbations for this 
system in one dimension. In two dimensions, however, the Calculus of Variations perturbations describe 
only non-dissipative perturbations. 
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The physical perturbations relevant to this system are incident charged particles and thermal fluctuations. 
They are not infinitesimal and certainly are dissipative. It would be more appropriate to allow perturbations 
in the normal current for infinitesimal calculations, and any finite perturbations calculated for two dimensions 
must use the TDGL in order to be relevant. This is the heart of the debate by Kramer et al. about whether 
perturbations on superheating can describe the first introduction of vortices into a superconductor. Just 
using the Landau-Ginzburg Hamiltonian is not appropriate to what would seem an immediately relevant 
problem. 

Perturbations using the equations of the Calculus of Variations approach should at least be close to 
the correct results. It is possible that the least stable subspace of dissipative perturbations is the space of 
stable perturbations. More ardent numerical work is required to test the dissipative perturbation equations 
which are of higher order than the conservative ones. We proceed with the conservative Landau-Ginzburg 
perturbations and their boundary conditions. 

The boundary conditions on the differential equations for the perturbations have to leave the applied 



fields unchanged. Given GL solutions like those shown in Figure 2.1 , the applied field, h(o) = q'{0) = HappUcd 
and /'(O) ~ must not change. Since the perturbation equations are linear, only the relative magnitude of 
the solutions is important. 

The eigenvalue, E, however, is nonlinearly coupled into the equation. However we scale the perturbations, 
the size of the eigenvalue will remain the same. It will depend only on the parameters (k, Ha) of the GL 
equation. In fact, the magnitude of the stability eigenvalue is another of the new results of this work. 
Knowledge of the eigenvalues as a function of perturbation wavelengths are useful in calculating higher- 
order stability. We begin as before with matched asymptotic analysis in inner and outer regions. 

2.3.2 One-dimensional Perturbations 

Outer solution. The outer equations for (/, q) are rescaled with before to yield (we will drop the 

subscript n for notational convenience) 

-/" + (3/2 + q'- 1)/ + 2fqq = Ef, (2.118) 



K^" + fq + 2fqf = Eq. 



(2.119) 



Expanding /, q, and E in powers of k, and recalling that q = to all orders in k in the outer region, we 
have at leading order 



-f'o + (3/o - l)/o = Eofo, 



(2.120) 



where /o — tanh . By changing variables to y — tanh ^ ^ ^ we see that the solution of Eq. ( 2.120| ) 

is the associated Legendre function of the first kind: 



foix') = coP^ 



tanh 



x' + xo 



(2.121) 



where /i = —^2(2 — Eq) and cq is a constant. The leading order solution for q is qo ~ 0. 

Inner solution. To obtain the inner equations, we rescale as in Eq. (2.66), with the perturbations rescaled 

as 



/V) = Fix), q{x') = K-^'^Q{X), 



(2.122) 
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such that 



\f" + (3^2 ^ 1q2 _ ^^p ^ = EF, 



K K 

?2 



-g" + F^Q + 2FQF = EQ. 



(2.123) 
(2.124) 



To leading order, Fq = 0, so that Fq = ao, with ao a constant. The leading order equation for the variation 
in Q is 



-Q2 + 2FoQoi^o + (i^o - ^^o)Qo = 0. 
The solution which satisfies the boundary condition Q'(0) = is 



Qo(^) 



2aoAo-Bo 



-AoX 



At 0(k) we find 



^1 = Qo-^o + 



(2.125) 



(2.126) 



(2.127) 



with the solution 

F^{X) = ai + aoBl 



Eo + AAl 



AAlEo 



-(Ao+s/ Al-Ea)X 



Eo {Ao + ^A'o-Eo)W^o-Eo 

4Al/Eo 



Eo + iAl 

2A0E0 {Ao + ^Al - Eo) V^g - Eo 



(2.128) 



X. 



We now have enough terms in the inner and outer region for a nontrivial match. 

Matching. We complete the matching of the inner and outer perturbations to obtain the eigenvalue, Eo- 
Performing a two term inner expansion of the one term outer solution, we have 



/o(«;X) ~ Co 



P!i{Ao) + ^sed,\xol^2)'^-^^^KX 
V2 



(2.129) 



where we have used tanh(a;o/-\/2) = ^o- Next, the one term outer expansion of the two term inner solution 



is 



Foix' Ik) + kFi{x' /k) ~ + 



ao2^'\l-Al) 



Eo 



Eo + ^Al 



AAl 



2 Ao (^0 + - Eo) ^/Ai - Eo 



(2.130) 



where we have used Bo = — 2^/''(l — Aq)^^^. Matching the two expansions using the van Dyke matching 
principle yields 



Co 



ao 



P2^{Aor 



(2.131) 



1 dP^{Ao) _ 2 
P^{Ao) dAo ^ e'o 



Eo 



■AAl 



AAl 



2^0 {Ao + - Eo) V^o - Eo 



(2.132) 



The last equation is a rather complicated implicit equation for the eigenvalue Eo{Ao), which generally must 
be solved numerically. However, when Ao = 1/V2 we find Eo = 0, corresponding to the critical case, with 
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-0.2 - 



0.0 0.2 0.4 0.6 0.8 1.0 



Figure 2.10: The stability eigenvalue E(Ao), with Aq the value of the order parameter at the surface 
at leading order. We see that E > for Aq > l/\/2, indicating locally stable solutions. 



> for ylo > 1/V2- The numerical evaluation of Eq. ( 2.132 ) is shown in Fig. 2.10 . Therefore, we see 
that our maximum superheating field (at lo west order) corresponds to the limit of metastability for these 
one-dimensional perturbati ons. In Fig. 2.11 we show Aq as a function of the lowest order magnetic field at 
the surface, Hq, from Eq. ( 2.9Cl| ). The stability analysis of this section shows that only the upper branch of 
this double valued function corresponds to solutions which are locally stable, with the field at the "nose" 
being the superheating field. 



2.3.3 Two-dimensional perturbations 

We next turn to the stability of the solutions with respect to two dimensional perturbations. It is very 
likely that there may be solutions stable with respect to one-dimensional perturbations but not two. The 
GL solutions are minimizers of the free energy and we expect them to usually sit in the well of free energy 
potential. They will likely always be stable with respect to infinitesimal one-dimensional perturbations. 
However, we can imagine that if we allow the free energy schematic a second direction, the GL free energy 
niinimizer may be either the minimum or maximum of a parabola in the y direction. In other words, we are 
searching for the applied field at which the GL solution becomes a saddle point in the free energy. 

If we perturb the extremal solution (/, q) of the GL equations by allowing f ^ f + 5f and q — > q + Sq, 
then the second variation of the free energy functional is 



S^J^ = dxdy 



^(V<5/)2 + 4/(<5/)q . <5q + fiSqf + {Zf + q^ - l){6ff + (V X Sc^f 



(2.133) 



We neglect perturbations along the i-direction because Fink and Presson showed that terms in z are 
purely positive definite and thus any variation in z only increases the free energy. Expanding in Fourier 
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Figure 2.11: The order parameter at the surface, Aq, as a function of the field at the surface, Hq, at 
leading order. The stability analysis shows that only the upper branch corresponds to locally stable 
solutions. The field at the "nose" is the limit of stability, and corresponds to the superheating field 
Ho = 2-3/4 ^ 0.595. 
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modes with respect to y , 

Sf{x,y) = f{x) cos ky, dq^{x,y) ^ qx{x)smky, Sqy{x,y) ^ qy{x) cos ky, (2.134) 

(0, (/(a;), 0), and integrating over y, we obtain (up to a 



substituting into Eq. (2.133), recalling that q 
multiplicative constant) 



6^ J' 



dx 



1 



1 



-^f ' + {if + q' + - l)f + Afqfqy + f\qi + q^) + {q'y - kq^Y 



(2.135) 



By integrating by parts and using the boundary conditions, Eq. (2.97), we can cast this functional into the 
form 



dxif,qy,qx)L2 




(2.136) 



where the self-adjoint linear operator L2 is given by 




1 



+ q^+ 3/2 + P/k' - 1 





(2.137) 



That the operator is self-adjoint shows that it represents perturbations on a conservative hamiltonian. As 
in the previous section, we want to determine the eigenvalue spectrum of this operator. We are primarily 
interested in the effects of long- wavelength perturbations (i.e., k 0), so we rescale k as k — nk' . Then the 
eigenvalue equations in terms of the outer coordinate (dropping the prime on k from now on) 



+ (3/2 + q^-l + e)f + 2fqq - Ef, 

-i^^Qy + f% + '^fQ.f - i^^kq^ = Eqy, 
n^kq'y + {f + ^k^)q,^Eq,. 

By using the last equation we may eliminate qx from Eq. ( 2.139| ), which becomes 



dx 



P-E 



p + K^k^-E^y 



+ fqy + 2fqf = Eqy 



(2.138) 
(2.139) 
(2.140) 



(2.141) 



For k = E qs. ( 2.1381) a nd ( 2.141 ) reduce to the one-dimensional perturbation equations of the last section, 
Eqs. ( 2.118 ) and ( 2.119 ); for E — the y redu ce to the Euler-Lagrange equations derived by Kramer [ p6| . 

The perturbation equations ( 2.138 ) and ( 2.141 ) may be solved by the method of matched asymptotic 
expansions, just as in the one-dimensional case. The derivation of the eigenvalue condition is essentially 
identical, with the final result that 



1 



dP^^iAo) 



P^{Ao) dAo 
V2(2 + So 



2 

E^ 



En 



AAl 



AAl 



2Ao 



[A^ + ^Ai^Eo)^A'^-Eo 



(2.142) 



where now /i — — ■\/2(2 + E^ — k^). The eigenvalue Ea{k) is plotted in Fig. |2.12 for several different values 
of Aq. For Aq > I/V2, Eo{k) > for all k, while for Aq < l/V^ there exists a band of long- wavelength 
perturbations for which Eo{k) < 0. In all cases the most unstable modes are at fc = 0, i.e., the one- 
dimensional perturbations are the least stable. This is in contrast to the large-K limit, where the most 
unstable mode occurs for fc 7^ G2I BO, El| . 



Landau-Ginzburg Systems February 1, 2008, 10:21 




Figure 2.12: The stability eigenvalue E(k) for two-dimensional perturbations of wavenumber k, for 
several different values of Ao. For Ao > l/\/2 the eigenvalue is stable for all wavenumbers, while for 
Aq < l/\/2 there exists a band of wavenumbers for which the solution is unstable. 



Landau- Ginzburg Systems February 1, 2008, 10:21 



38 



10" 



10' 



10 



lo-" 



10' 




Figure 2.13: This figure shows the numerically calculated superheating field for large ft (solid line) 
compared with the two-term asyptotic expansion derived by Chapman (dashed line). The slope of the 
dashed line is —4/3. 



2.3.4 Large-K two-dimensional Stability 

The exact numerical solution of the Landau-Ginzburg equations shows that the system is unstable with 
respect to two-dimensional perturbations at an applied field, H2B < Hsh. Fink and Presson |Q estimate 
that the H2D separates from Hgh at k « 1.10 or k « 1.13. Our calculations show the crossover occurs at 
1.16 < K < 1.17. At the bifurcation, the least stable mode lifts from A: = to steadily higher wavenumbers. 

Chapman pl| ] recently used matched asymptotics to examine the Landau-Ginzburg equations in the 
high-K limit. His final result for the superheating field is 



Hsh — 



1 

V2 



(2.143) 



where the constant C is determined from the solution of the second Painleve transcendent; a numerical 
evaluation yields C=0.326 The first term was originally derived by Ginzburg Q, and the second term 
with the unusual dependence on k is new. We verify the second term in Fig. 2.13. 
but does not converge as rapidly as the small-K solution. 



The dependence is correct 



2.4 Numerical methods 



We used two separate algorithms to evaluate our solutions of the superheating field. First, we calculated 
solutions to the Landau-Ginzburg equations for the half-space as a function of (k, iJappiied)- Second, we 
calculated perturbations to those solutions as a function of wavenumber, fc. 

Finding solutions to the Landau-Ginzburg equations themselves was straightforward. To ensure that no 
current passes through the boundary at a; = and that the sample is totally superconducting infinitely far 
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from the surface, we impose the boundary conditions 

/'(O) = 0, fix) -> 1 as X ^ oo. (2.144) 

Since the field at the surface must equal the applied field Ha, and there must be no field infinitely far from 
the surface, we impose the boundary conditions 

h{0) = Ha, q{x)-^0 asx-^oo. (2.145) 

The discretization used a finite domain, so the boundary conditions at infinity were generally enforced at 
a coordinate large enough not to change the shape of the solutions of the independent variables (/, q) as 
measured by the norm of the difference between successive solutions, 

J i\fl{x) - /2(x)| + \qi{x) + q2{x)\) dx. (2.146) 

Later investigations instead used boundary conditions which were the analytically derived asymptotics of 
the Landau-Ginzburg equations. They had no measurable affect on the independent variables. 

For K — !■ 0, we rescale the equations as x' = kx making the new unit of length the correlation length ^. 
Since ^ ^ A in this limit, a numerical solution over a domain much larger than ^ would ensure that the 
regions of rapid change for / and h would be included. (For small k, we find that solving for x' < 500 is 
sufficient.) In the large k limit, we use the rescaled equations again, but we increase the size of the domain 
depending on the value of k. (The equations must be solved for domains as large as x' < 10^ for values of 
K ~ 103.) 

The equations can be solved using the relaxation method By replacing these ordinary differential 
equations with finite difference equations, one can start with a guess to the solution and iterate using a 
multi-dimensional Newton's method until it relaxes to the true solution. In order to more accurately pick 
up the detail near the boundary, we choose a grid of discrete points with a higher density near x — 0. In 
particular we choose a density which roughly varies as the inverse of the distance from the boundary. (For 
low K our density, in units of mesh points per coherence length, varies approximately from 10^ near the 
boundary to 10^ at the farthest point from the boundary, while for high k it varies from 10^ to 10^^.) 

Hsh can be found in the following way. For a given value of k an initial guess is made where there 
is no applied field and the sample is completely superconducting {f = l,q = 0,h = 0). The field Ha is 
then increased in small increments. For each value of Ha a solution is sought using the result from the 
previous lower field solution as an initial guess. Eventually a maximum value for Ha is reached, above 
which one of two things happens: our algorithm fails to converge to a solution or it converges to the normal 
(nonsuperconducting solution). This maximum value of Ha is the numerical result for H^h- Using this 
algorithm, Hshin) can be found for a wide range of k's. 

It is possible to imagine a situation in which this algorithm might not work. Suppose that you have 
a solution at Hai and fail to converge on a solution at Ha2 ~ Hai + 0.01. It is possible that your initial 
guess was just not close enough. For instance, a smaller stepsize would permit one to find a solution first 
at Hal + 0.05, and then that solution would be close enough to find the solution at Ha2 — Hai + 0.01. 
The superheating solution could creep away as quickly as you could approach it. Using a variation on the 
Contracting Mapping Theorem, Herbert Keller |^ showed that Euler's Method, in particular, guarantees 
that there exists a finite neighborhood of the solution which will always converge for well-behaved systems 
such as ours. 

Each run (for a given n) takes about 10 seconds on a Pentium II 300. We find it sufficient to deal with 
superheating field values for 10"'^ < k < 10^. 

More interesting is the study of the perturbations. These require a third parameter, the wavenumber 
k. The solutions are in a three-parameter space {K,Ha,k) and the solution of each perturbation requires 
calculation of the initial Landau-Ginzburg solution. Because the two-dimensional perturbations were of 
interest for large n where there are not readily available analytic solutions, even asymptotics, the initial 
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conditions for each run (of wavenumbers) are finicky at best. 

The primary objective of solving the two-dimensional perturbations is to find the H2B line, but solutions 
to our equations yield, in addition, the dependence of the stability eigenvalue on the wavenumber. That is 
more information than is relevant to the superheating field, so it is not included here. 

Calculations of the two-dimensional stability required significantly more computer resources, both in 
CPU time and storage space. The algorithms were written to function on a cluster of a dozen systems and 
save only relevant data in indexed, binary files on a central server. Time-critical sections were re- written as 
FortranQO subroutines to C"*""*" control structures. Finding the bifurcation point at k ~ 1.17 required about 
a week on our Pentium II cluster or three months of computer time on a single machine. 

2.5 Discussion 

We have solved the Landau-Ginzburg equations both analytically and numerically for a superconducting 
half-space. The asymptotic methods depend on disparity between coherence length and penetration depth, 
but the solutions remain relevant even where they are equal. The resulting expansions for the superheating 
field should be immediately useful for the reverse operation — calculating the Landau-Ginzburg parameter 
from the superheating field. 

The same techniques were also effective for deriving perturbations on the superheating field. These were 
previously deemed complicated enough to be beyond analysis by most authors. We not only calculated the 
exact two-dimensional perturbations but also elucidated a vexing question last posed by Kramer pof about 
whether perturbation solutions can represent vortex nucleation. 



Chapter 3 



Phase Transition in a 
Current-carrying Wire 



3.1 Introduction 

When a superconductor is placed in a magnetic field equal to its critical field He, the normal and supercon- 
ducting phases can coexist in a state of equilibrium with the two phases separated by normal-superconducting 
(NS) interfaces. The dynamics of such interfaces is important for various nonequilibrium phenomena. For 
instance, if the applied magnetic field is quenched below He, these interfaces move through the sample, 
expelling the magnetic flux so as to establish the Meissner phase |7^, ^ 18, ^ |l^, |5^. Just as 



superconductivity can be destroyed by applying a magnetic field exceeding He, it can also be destroyed by 
applying a current exceeding the critical depairing current Jc- Thus by analogy with the magnetic field case, 
one might expect the competition between the superconducting sample and the applied field to stabilize an 
NS interface in a current-driven system [Q. In contrast to the magnetic field induced NS interfaces, these 
current-induced NS interfaces are intrinsically nonequilibrium entities, and their structure depends upon the 
dynamics of the order parameter and magnetic field. The evolution and dynamics of these nonequilibrium 
interfaces is the subject of this chapter. 

The current-induced NS interfaces arise in several contexts. First, they are known to be important in 
understanding the dynamics of the "resistive state" in superconducting wires and films (for a review see 
p7| or |lOl| ), and in determining the global stability of the normal and superconducting phases in the 
presence of a current [|67| . Second, Aranson et al. [0 have recently used simulations to study the nucleation 
of the normal phase in thin type-II superconducting strips in the presence of both a magnetic field and a 
transport current. They found that a sufficient current produced large normal droplets containing multiple 
flux quanta. Without a current one finds stationary, singly quantized vortices, with a larger amount of NS 
interface per flux quantum than a multiply quantized droplet. They conclude that the current produces an 
effective surface tension for the NS interface which is positive, stabilizing the interface and producing larger 
droplets with smaller surface area. In their simulations, topological singularities of multiple flux lasted for 
the duration of simulations of flux entry from demagnetizing fields. Motivated in part by its role in this 
phenomenon we wanted to re-examine the nonequilibrium stabilizing effects of current. 

Even when the superconducting phase is ostensibly the equilibrium phase, a current makes the normal 
phase metastable, i.e., linearly stable to infinitesimal superconducting perturbations. A localized supercon- 
ducting perturbation of finite amplitude, on the other hand, has one of two fates: (1) its amplitude may 
ultimately shrink to zero restoring the normal phase (undercritical) or (2) it may grow eventually estab- 
lishing the superconducting state (overcritical). Separating these two possibilities are the critical nuclei or 
threshold perturbations, for present purposes stationary solutions of the time-dependent Ginzburg-Landau 
(TDGL) equations localized around the normal state. As one raises the current, the amplitude of the 
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normal superconducting 

~1 

globally stable metastable 

J* 

metastable globally stable 
J = 

Table 3.1: We examine transitions between the two homogeneous states of a current-carrying super- 
conductor. While the system docs not conserve energy, we can define metastability with respect to 
small perturbations. 



threshold solution grows, implying that the normal phase becomes increasingly stable. 

At very low currents, the widths of the critical nuclei shrink as the current is increased, but eventually this 
trend is reversed and the width grows as the current is increased further. In fact, as the current approaches 
a particular value, J* (the "stall current" the width diverges resulting in two well-separated, stationary 
NS interfaces. Above J*, no nucleus solutions exist in the TDGL. The absence of a critical nucleus for the 
superconducting state defines the normal state above J* to be globally stable |6^, |3^ . 

When a thermal fluctuation is larger than the critical nucleus for a particular critical current, that 
perturbation will grow first to locally saturate the order parameter, then form separated NS interfaces 
which move at constant velocities controlled by the applied current. The transformation of random thermal 
fluctuations into NS interfaces of fixed form is phase ordering. Once an interface forms, it will travel towards 
the normal phase if J < J* and towards the superconducting phase if J > J* . 

The interface solutions have been studied numerically by Likharev |Q, who found that the interfaces 
were stationary at J* w 0.335 for u — 5.79, where u characterizes the material and is 5.79 for nonmagnetic 
impurities |Q. They were also studied by Kramer and Baratoff |6^, who found J* « 0.291 for u — 12 
(corresponding to paramagnetic impurities pl|]). However, we know of no systematic study of the dependence 
of J* upon u. In this work we remedy this situation by using a combination of numerical methods and analysis 
including matched asymptotic expansions fl^, 105 . We show that J* ^ u~^^^ for large u in contrast to a 



previous conjecture |74[] , and we find how J* approaches Jc in the small-u limit. 

At currents close to J* , we can treat the interface velocity as proportional to {J— J*) and calculate a kind 
of susceptibility. Likharev [f74| defined the constant of proportionality, rj = {dc/ dJ)~^\j^j- , where c is the 
interface speed; he found rj w 0.7 for u =^ 5.79. In the extreme limits, J — > and J Jc Likharev predicted 
that the speed c diverges. We find c to be bounded in both cases and provide an analytic expression for it 
as J ^ 0. 

The results of this work are summarized in Table 3.2. The rest of the chapter is organized as fol lows. 
After briefly reviewing the TDGL equations and the approximations used in this work (section ^.3.1 ), we 
study the critical nuclei focusing on their size and shape in the limit J — > (section |3.4D . We then move on 
to consider the stationary interface solutions; in particular we map out the dependence of the stall current 
J* on u and supplement the numerical work with analysis of the u — s- oo and u limits (section ^.5[ ). 
Next, we examine moving interfaces first in the linear response regime and then in the limits J — > and 
J ^ Jc (section 3.6). Appendix |^ contains a calculation of the amplitudes of the critical nuclei in the J ^ 
limit. 



3.2 The Physical System 

Our simplified picture of a one-dimensional superconductor carrying a current clarifies some of the essential 
physics of more complex phenomena in experiments on thin whiskers or strips of superconducting mate- 
rial. A general review of current-induced phenomena in one-dimensional superconductors can be found in 



Tidecks [101|. In this section, we will discuss properties of samples (dimensions, materials, contacts) and 



give a brief taxonomy of behaviors. 
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Table 3.2: Summary of the primary results. 



I. Critical nucleus 

Small- J width 
Small-J amplitude 

II. Stall current J* 
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Sec. 
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Figure 3.1: A phase-slip center is formed when a local perturbation of the current causes the order 
parameter to drop to zero. When the order parameter is small, the phase of the order parameter 
twists another 27r. The order parameter then heals slightly leaving a localized area with higher normal 
current and a supercurrent decreased to about 25 % of the total current. 



Systems are quasi-one-dimensional because the coherence length is larger than the lateral dimensions 
of the superconducting material. Both whiskers of single crystal and etched film depositions can meet this 
criterion easily. The point is that the material is thin and skinny enough that variations in the order 
parameter and its phase are not significant across the width of the superconductor. A sample of YBCO from 
Jelila et al. was 200 ^m long, 20 /^m wide, and 90 /ini thick. This precludes the formation of vortices in 
the sample by self-induced magnetization from applied currents. 

We have discussed already the metastability of the two uniform states, the normal and superconducting 
states. Measurements in steady state rarely show a transition from normal to superconducting or vice versa. 
Instead, the system passes into an intermediate state called the resistive state where the sample is mostly 
superconducting but there are occasional slips in the order parameter. These slips are oscillatory regions 
where the order parameter decreases and more of the current is carried by both an increase in t he c hange 
of phase of the order parameter (J, cx V'* VV^ — i/'V?/'*) and by a localized normal current. Figure 3T shows 
a single phase slip. At currents above J*, the resistive state consists of a periodic array of these oscillatory 
regions, called phase slip centers, PSC. The PSC state allows the strip to remain superconducting above its 
critical current, Jc- 

When the current rises above the critical current, a single PSC will enter. More PSC enter at successively 
higher applied currents. The result is a series of steps in the I-V curve called a "forked ascension," shown 
in Fig. |3.2| . Experiments on tin |98, Q and YBCO ||6l| demonstrate similar voltage characteristics. Each 
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Figure 3.2: These are the dc current- vohago characteristics of a YBCO bridge from Jelila et al. | |6l| . 
The dotted hues show hysteresis for rising and falling currents. This I-V is typical of pure supercon- 
ductors not too far from transition temperatures. 
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Figure 3.3: Even when the system switches from the normal to the resistive state, the process should 
occur by means of a moving phase front. 



step represents the addition of a PSC. The first sohd analysis of PSC by Skocpol, Beasley and Tinkham ||9q| 
attributed the stabihty of the centers to thermal effects. Because the middle of the phase slip is normal, 
it heats, possibly above the transition temperature. Averaging the additional resistance of each PSC, they 
used the normal state resistivity to find that the length scale of a PSC is the quasiparticle diffusion length. 
They mention, however, that in very clean systems near Tc, steps are still observed "due to some intrinsic 
inability of the uniform state to nucleate in the presence of the phase-slip process." There is debate over 
when the normal or superconducting systems are unstable with respect to the formation of phase slip centers, 
but they appear to be intrinsic to the isothermal behavior of the order parameter and a necessary state in 
the adiabatic transition from normal to superconducting. 

This does not mean, however, that a transition from normal directly to superconducting is not possible. 
There are two ways to see the transition directly. As is shown in in Pals and Wolter and Jelila et 
al. | |6^ , the order parameter can take up to 200 ns to respond to a 3 ns change in the applied current. It 
is possible to drop the current below the PSC nucleation current before the system can respond. Secondly, 
applied currents in a system near Tc can be so small a PSC cannot form. Each phase slip center has a kind 
of activation energy, the amount of excess current required to depress the order parameter enough that it 
reaches zero so a phase slip can occur. If the critical current is less than the amount of current required to 
form a PSC, then the transition will be direct from normal to superconducting. 

Also, if a sample has superconducting contacts, it will be more likely to become superconducting at higher 
temperatures and higher applied currents \ ^t\ . While most experiments etch bridges into uniform substrates, 
some have applied normal contacts to thin strips These would be more amenable to measuring direct 
phase transitions .[| 

The transition from normal to the PSC state is not so different from the transition from normal to 
superconducting. While, in decreasing currents, the PSC state forms much the same way the Abrikosov 
lattice forms in a type-ii superconductor, the system's reaction to a sudden decrease in applied current 
would be nucleation. The supercurrent at the contacts nucleates a local phase change which spreads across 
the superconductor switching it from normal to the resistive state. The initial wavefront is probably just 
a switching wave and is exactly what we study here. The subsequent relaxation to the resistive state is 
not examined here, however. A simple picture is shown in Figure [3.3| . In any case, there seems to be little 
experimental data on this region of the phase diagram. 

Even if one can avoid forming PSC, heating of the normal metal masks some of the intrinsic supercon- 
ducting properties of thin superconductors. Critical currents of cold superconductors can be large enough 
that resistive heating in normal regions can control the stability of stationary structures ||9^ and the veloc- 
ity of NS boundaries . Because we are interested in the coexistence of superconducting and dissipative 
normal regions, heat flux can be a significant factor in thermally isolated systems. 

There is a large literature devoted to heat flux in current-driven superconductors. Hot spots can form 
stable autosolitons or drive phase transitions for the entire system. We are more interested in the behavior 
of isothermal systems because they describe fundamental behavior of superconductivity. The velocity of a 
normal-superconducting interface, for instance, will be limited by some combination of heat flux radiation, 

^ Supercondu cting contacts measure only the pair chemical potential while normal contacts measure the full potential of the 
electric current [ 106 . If the edges are in steady state, however, the difference should be negligible. 
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diffusion of the order parameter, and relaxation times of superconducting pairs into quasiparticles. The 
experiments and analysis presented here will be related to systems shown to be highly isothermal so that 
the rate of heat flux will be of minimal importance. 

There two ways commonly mentioned to combat heating effects. In the second of their set of articles on 
nonuniform states in thin strips in 1974 |99[ , Skocpol et al. point out that near Tc critical currents are small 
so that applied voltages will heat the sample less. More often, experimenters hope that a superconducting 
strip on a substrate with good thermal capacitance will dissipate heat efficiently. 

More enlightening are time-resolved studies of voltage as an applied current pulse is modulated from 
above the critical current to currents below J* and Jmin- One experiment which seems a clear attempt to 
measure phase changes at applied currents below the critical current is Jelila et al. |6^ . The measurements 
are so appropriate they must be discussed despite inconclusive results. They applied a varying current pulse 
to a thin strip of YBa2C302, in order first to drive it normal then to reduce the current and watch it return 
to the superconducting or PSC states. While their main focus was an initial time delay of the material in 
responding to applied currents, this delay was much better explained in a later PRL by the same group 
The earlier paper shows two graphs of voltage versus time with resolution in nanoseconds. 

The film itself was grown on a MgO substrate which was verified to have very efficient heat conductiv- 
ity 1?^. The film was a 30 nm thick, 200 /im long, and 29 /zm wide microbridge. The normal state resistance 
of the bridge seems to be around 5 il. Experiments were conducted at 4.2 K. We discussed earlier that 
temperatures close to Tc are better suited to isothermal measurements, but the signal to noise ratio was 
limited by shunt resistors. 

The paper shows two experiments on slightly different samples. The first applies a supercritical 67 mA 
pulse which drops to 20 mA. Here, the voltage falls within the fall time of the pulse generator, 3 ns. If we try 
to explain the phase change as the progress of wave fronts, we can conclude the wave fronts traveled faster 
than 3.3 x 10* m/s. 



The second experiment is shown in Figure 3.4. While the initial current of 70 mA drives the system 
normal we see a broad inductive rise in the voltage. Then the applied current is lowered to less than the 
critical current, and the system shows one of two behaviors. If the applied current is greater than 57 mA, 
the system relaxes into a PSC state with a voltage of 92 mV. If the applied current is less than 57 mA, the 
system returns to the superconducting state. The authors say only that "the resistance drops to zero in 
about 65 ns." The shape of the curve is interesting. At early times, we expect the shape to be dominated by 
rapid phonon removal. Later times, however, show a curious shape of the curve which is not inductive. 

It is unclear whether the system has had enough time to form the PSC state before undergoing a transition 
into the pure superconducting state. Suppose that the last downturn of the voltage, from 207 ns to 224 ns, 
is related to the movement of normal-superconducting phase fronts in the material. The shape of that part 
is roughly linear. The change of resistance with time of 0.0742 fi/ns suggests phase front velocities of about 
1600 m/s which is in the general range of the calculations in this chapter. 

What exactly is occurring during the phase change is not at all clear. There doesn't seem to be much 
data on phase changes at nonzero currents. Rather, most experiments examine the rate of phase transition 
at zero bias or currents greater than the critical current. 

3.3 Models for Phase Transition in a Wire 
3.3.1 The TDGL Equations 

The starting point for our study is the set of TDGL equations for the order parameter tp, the scalar potential 
<&, and the vector potential A: 



a|V-6|^|V, (3.1) 
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Figure 3.4: This shows the voltage response as a function of time for a YBCO microbridge from 
Jelila et al. |^^. The higher voltage drives the system normal. The system then returns either to the 
superconducting state, Va, or the PSC state, V5. 



47r 

V X V X A = — (J„ + J,), 

c 



where the normal current J„ and the supercurrent are given by 



0*2 



he* e" 

— (V^* - ^W/) - — IVI'A, 

Zmi mc 



(3.2) 

(3.3) 
(3.4) 



and where 7 (which is assumed to be real) is a dimensionless quantity characterizing the relaxation time 
of the order parameter, cr(") is the normal conductivity, and a = ao(l — T/Tco). From these parameters 
we can form two important length scales, the coherence length ^ — ?i/(2m|a|)^/^ and the penetration depth 
A= [m6cV4^(e*)'|a|]i/2. 

These equations assume relaxational dynamics for the order parameter as well as a two-fluid description 
for the current. With somewhat restrictive assumptions, they can be derived from the microscopic BCS 
theory ]9^ , Further simplification is possible in the limit of a thin, narrow film, that is, when the 

thickness is less than the coherence length, d < ^, and the width is less than the effective penetration depth 
p9[ , w <^ X^/d. In this case the current carried by the film or wire is small, and we needn't worry about 
the fields it produces. Therefore, Eq. (3.2) may be dropped, and we need only specify the total current 
J = J„ + Js (subject to V • J = 0), along with the order-parameter dynamics, Eq. (3.1). This approximation 
is commonly used for superconducting wires p7[ | and can be justified mathematically for superconducting 
films pO| . In addition, we will be considering processes in the absence of an applied magnetic field, so that 
we may set A = 0. With these simplifications, we can now rewrite the equations in terms of dimensionless 
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(primed) quantities, 



e*^|a| 

^=V^^A (3.5) 



which leads to 



u{dt' + = (V'2 + 1 - )V'', (3.6) 

J' = Im(V/*V'V/) - Vy, V'-J'==0. (3.7) 

Note that length is measured in units of coherence lengthQ. We will drop the primes hereafter. The only 
parameters remaining in the problem are the scaled current J and a dimensionless material parameter 
u = Ttp/rj, where = ?i7/|a| is the order-parameter relaxation time and tj = a^-^^^mb/e*'^\a\ is the current 
relaxation time. We will treat u as a phenomenological parameter and study the nucleation and growth 
process as a function of u. The microscopic derivations of the TDGL equations predict that u — 5.79 
(nonmagnetic impurities) ]96| , and u — 12 (paramagnetic impurities) but small u is also useful for 
modeling gapped superconductors | |60[ |. 

3.3.2 Generalized TDGL Equations 

There have been fruitful generalizations of the TDGL applied to superconducting wires. The simplest 
recognized that the change in the magnitude of the order parameter is more closely related to the relaxation 
time of the superconducting pairs while the change in the phase of the order parameter is more closely related 
to relaxations of the quasiparticle excitations. These authors substituted two relaxation constants, and 

More fruitful were several competing derivations of TDGL which account for superconductivity with 
a finite gap. The first was Kramer and Watts-Tobin Bq], but alternative derivations and discussions are 



available in Ivlev and Kopnin (57| and Tidecks [101|. All derived from microscopic theory, each accounting 
for somewhat different pair-breaking mechanisms which lead to different contributions to a gap parameter, 
r. The general form of the equations is 

-u + 1) + + 2r^) + W + - m = (3.8) 

j = - V0 + 1 (V/ V^j - ^pvr )■ (3.9) 

When the gap parameter is zero, these equations reduce to the standard TDGL. 

Having more parameters clearly promises greater specificity to particular metals, but the gap parameter 
explains significant shortcomings in the how the TDGL describe the PSC state. Whereas we show below 
that the basic TDGL do not allow a normal to superconducting transition above J*, the equations with a 
gap do demonstrate such a transition. They also better explain the experimentally observed stability ranges 
of the resistive state. 

The most creative and appropriate generalization of the TDGL for this system is in a paper by Eckern, 
Schmid, Schmutz, and Schon [|33| . They derive a nonlinear Langevin equation similar to the TDGL which 
is appropriate to dissipative superconducting systems. More interesting is that they are able to develop a 

^We warn the reader that this choice is different than in many applications of the TDGL equations, where the penetration 
depth is chosen as the length scale. In the thin film limit the penetration depth drops out of the calculation, leaving ^ as the 
length scale. 
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measure on this system that aUows them to work with the Langevin equation much hke one treats a typical 
conservative free energy. They then proceed to describe the basic system stabihty in terms of metastable 
and stable configurations and calculate threshold solutions from critical nuclei, much as is done here. It is a 
fascinating attempt to coax sensible physics out of a nonlinear dissipative system. 

3.3.3 Heat Equations 

As mentioned earlier, there is a large literature devoted to the analysis of normal-superconducting boundaries 
driven by heat flux. This analysis is entirely appropriate to a large class of thermally isolated systems. It 
seems most appropriate, for instance, to current flow through a wire. The analysis of these systems also 
leads to coupled nonlinear equations, usually parabolic diffusion equations. They also display traveling 
autosolitons as well as a wealth of other solitons common to reaction-diffusion equations. 

The first work on thermal effects in domain boundaries was done by Skocpol, Beasley, and Tinkham ]9^ , 
[ggf . The did steady-state calculations for heat flow in thin bridges where they assumed the normal- 
superconducting boundary was sharp. They predicted stable solitons where the center of the bridge remained 
normal while the edges were superconducting. 

Recently, Rudyi has both re-examined the stationary thermal solitons and switching waves [ p4| , 
which we study here for isothermal systems. When looking at thermal variations, one uses the temperature 
as an independent variable rather than the order parameter. Using the notation, Q{9) = {T{9) — To)/Tc, 
where T{9) is the film temperature, 9 — x — vt a, self-similar variable, Tq the coolant temperature, and Tc 
the critical temperature, Rudyi models the system as 

Qm + -Q's{0)~h,Qs{9)^G (3.10) 

e'^(0) + —Q',M - bnQn[9) + W^0. (3.11) 

The constants concern heat transfer and thermal conductivity. These equations predict a bistable, hysteretic 
system which collapses through critical nuclei which become phase fronts. While thermal propagation is 
a different mechanism from isothermal phase propagation, the behavior of the system is almost identical. 
Most thermal models are well understood as classic reaction-diffusion systems |55| and, as such, they can be 



analyzed with the standard autosoliton theory [103, 104| . As shown in Appendix the TDGL do not quite 
conform. 

While thermal flux can constrain the velocity of a normal-superconducting interface or create nuclei in a 
filament, it's effects can be minimized by good thermal coupling of the superconductor to a substrate. The 
order parameter relaxation and current relaxation are then the dominant factors in the stabilization of the 
normal-superconducting interface or the growth of nuclei. 



3.4 Nucleation of the Superconducting Phase from the Normal 
Phase 

In the presence of an applied current the normal phase in a wire is linearly stable with respect to supercon- 
ducting perturbations for any value of the current [t^ . The reason for this stability is that any quiescent 
superconducting fluctuation will be accelerated by the electric field, its velocity eventually exceeding the 
critical depairing velocity, resulting in the decay of the fluctuation. The growth of the superconducting 
phase therefore requires a nucleus of sufficient size that will locally screen the electric field and allow the 
superconducting phase to continue growing; smaller nuclei will simply decay back to the normal phase. The 
amplitude of the "critical" nucleus should decrease as the current approaches zero, reaching zero only at 
J — 0. We expect the critical nuclei to be unstable, stationary (but noncquilibrium) solutions of the TDGL 
equations, which asymptotically approach the normal solution as x ±oo. These "bump" solutions of 
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the TDGL equations are the subject of this section. We include here an extensive numerical study of the 
amplitudes and widths of the critical nuclei, as well as some analytical estimates for these quantities. 

3.4.1 Numerical Results 

Let us start by discussing the numerical work on the critical nuclei. For the analytic work, we often find it 
convenient to use the amplitude and phase variables, i.e. tp = /e*^; but they are ill-suited for the numerical 
work, since the calculation of the phase becomes difficult when the amplitude is small. Following Likharev 



JtJI we use instead = R + il , with R and / real, and in one dimension Eqs. (3.3.1) become 

uRt ^ R^^+ufil + R- {R^ + I^)R, (3.12) 
ult = I^^-ufiR + I~{R^+I^)I, (3.13) 
J = RI^-IR,-fi^. (3.14) 

Since the nuclei are unstable stationary states, they are investigated only by time-independent means. 
Such solutions require a particular gauge choice — in this case fi{x) = where ip{x) has its maximum ampli- 
tude; they are then sought using a relaxation algorithm p3]. Figure shows a typical bump's amplitude. 



/ = VRF+I^, the associated superfluid velocity q = {RIx — IRx)l P and the electric field E(x) = —iix{x). 
The figure shows only half of the solution; /(x), q{x) and E{x) are even about a; = 0. In Fig. |3.6| we plot the 
bump's maximum amplitude, ipo, as a function of J; it grows as the current rises, indicating the increasing 



stability of the normal phase. In the data presented by Watts- Tobin et al. [106| V'o appears to vary linearly 
with J for small J. However, in our numerical calculations at very small currents the dependence deviates 
from linearity (see the inset of Fig. 3.6|) , and V'o drops rapidly to zero as J ^ 0, consistent with the exponen- 



tial behavior suggested in Refs. ||57, |59|. More precisely our small- J data (0.008 < J < 0.015) at u = 5.79 is 
fit by 

ipoi-J) ^ Bexp{~A/uJ), (3.15) 

with A — 0.042 and B = 0.19. A somewhat similar dependence (with A — 2/3) was suggested by Ivlev et 
al. Il^, |5^ ; they were considering a distinct quantity but one also related to critical fluctuations about the 
normal phase (see the Appendix for more details). 

The width of the bump diverges in the small-,/ limit like (mJ)~^/^, as can be seen from the analysis 
below. So as J is increased from zero, the width initially shrinks, but eventually the width begins to grow 
again, diverging as the current approaches the stall current J* . In this limit the bump transforms into two 



well separated interfaces (see Fig. 3.7) 



3.4.2 Analysis in the J — Limit 

The equations for nuclei centered at the origin are 

^/'^^-mW-' + V'-|V'PV' = 0, (3.16) 
li, = —Jx+ I Im ('0*-0a;/) dx', (3-17) 







where we have dropped the term ipt and selected the gauge /i(0) = 0. We saw in Fig. 3.6 that ■00 becomes 
very tiny in the small-J limit, thus the nonlinear terms can be neglected, leading to 

%pxx + iuJxTp + ip — 0, (3.18) 
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Figure 3.5: The bump's amplitude, f{x), its superfiuid velocity, q(x), and the electric field, E{x), for 
u = 5.79 and J = 0.2. 
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Figure 3.6: The maximum amplitude of the bumps as a fu nctio n of J for u = 5.79. The inset 
shows the exponential dependence of the small-J data, see Eq. (B.lSh. 
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a complex version of the Airy equation. Applying the WKB method results in the approximate solution 



[1 + [uJxf]-^/^ 
2 



exp 
exp 



3 uJ 



21 3/4 



3uJ 



[l + (uJx) 



a 



(3.19) 



where a = tan ^{uJx). The numerical data agrees quite well with this predicted shape in the small- J limit. 
For small x the expression can be approximated by 



tp ~ exp - uJ/A)x - uJx"^ /A\ 



(3.20) 



We see here that the width of the bump varies like {uJ) in this limit and that the superfluid velocity 
g « (1 — uJ/4). For large x, on the other hand, where a « 7r/2, the expression becomes 



^ {uJx) exp 



V2itJ, 



:|3/2n - 



(1-i) 



(3.21) 



as one expects for the Airy function. Note that deep in the tail of the solution, we see a different length 
scale \Airy ^ {uJ)~^^^ arising. 

Since the; above analysis is of a linear equation, it can not determine the amplitude of the nucleus; for 
this purpose the nonlinearities must be considered. In the appendix we outline an ad hoc calculation of the 
small- J limit of the bump amplitude. We take a of an unknown amplitude but of a fixed shape inspired 
by the above analysis and assume that it is a stationary solution of the full TDGL. We then determine its 
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amplitude sell-consistently. The resulting amplitude is 



2J\'/V9 1\"'/' / 16 



1^1 "UJ [s-u) ^^n^SnijJ- ^'-''^ 
The lactor A — 16/81 is within a few percent of that extracted from the numerical data. 

3.5 Stationary Interfaces 

As the current is raised, the width of the critical nucleus grows and ultimately diverges as the stall current 
is reached, resulting in well separated, stationary interfaces. These interface solutions will be the subject of 
the rest of this work. 

3.5.1 Numerical Methods and Results 

Let us first discuss the numerical work on the interface solutions. For given values of u and J we evolved 
the TDGL equations from an initial guess which is purely superconducting on the left, ip{x) — foo e"^°°^ and 
l^x{x) = 0, and purely normal on the right, •4>{x) = and ^J-xix) = —J. The values foo and goo are related 
to the applied current through 

J = f lV^~ fE: (3.23) 
qoo = yr^. (3.24) 




Stability requires taking the larger positive root of the former equation 1 73 which places the following bounds 
on J, foo and goo: 

v/4/27 w 0.3849, (3.25) 
2/3 « 0.8165, (3.26) 
T73w 0.5774. (3.27) 

We employed several schemes to integrate the equations in time including both explicit (Euler) and implicit 
(Crank-Nicholson) |9|. 

Initially the front moves and changes shape but eventually it reaches a steady state in which the interface 
moves at a constant velocity without further deformation. By the time-dependent means we found locally 
stable, constant-velocity solutions for currents less than Jc- To examine these solutions more accurately, we 
adopted a time-independent method. First, we transformed coordinates to a moving frame, x' — x — ct] 
next, we chose /i = cgoo as a: ^ —oo which allows for a truly time-independent solution (i.e. one with 
both amplitude and phase time- independent). Then we searched for stationary solutions using a relaxation 
algorithm where (u, J) are input parameters and c is treated as an eigenvalue. This approach requires 
an additional boundary condition to fix translational invariance; we elected to fix /i on the rightmost site. 
To find the stall currents J* we can set c = and take J or u as the eigenvalue. 



Figures 3.8 and 3.9 show the order-parameter amplitude / and the electric-field distribution E = —^Xx of 
the stalled interface determined for u — 500 and u — 1.04 respectively. Note that while / is very flat in the 
superconducting region, the real and imaginary parts, R and /, oscillate with a wavelength 27r/goo- Because 
of this additional length scale inherent in R and /, there is little to be gained from varying the mesh size. 
In fact, this length scale is compressed as we move to the right, and we are only saved from the difficulties 
of handling rapidly oscillating functions by the fact that the amplitudes decay so quickly. 



In the large-u case (see Fig. 3.8), E(x) remains flat throughout most of the space; it changes abruptly from 
one constant to another only after f{x) has become small. The variations of f{x) are more gradual; however, 
the greatest changes in fx occur in that same small area. This region of rapid change is known as a boundary 
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X 



Figure 3.8: The stationary NS interface solution when u = 500 for which the stall current J* = 0.12252. 
Shown here ar e the numerically determined f{x) and E{x), as well as the Langer-Ambegaokar (LA) 
solution (Eq. (3.43), the solution with no electric field) corresponding to the same current. 




X 



Figure 3.9: The stationary NS interface solution when u = 1.04 for which the stall current J* = 0.3836. 
Shown here are the numerically determined f{x) and E{x), as well as the function /o(i')> t^*^ « — + 
profile, derived from Eq. ( 3.7C ) where x = u^^^x. 
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layer, it marks where the current suddenly changes from superconducting to normal, i.e., the position of the 
NS interface. As u increases, the longer length scale over which / varies on the superconducting side remains 
essentially fixed, while the boundary-layer thickness shrinks to zero. In the opposite limit, the small-u case 



(see Fig. 3.9), f{x) and E{x) appear to vary together even in the superconducting region; moreover, the 
length scale over which they vary grows as u is decreased. We will postpone providing more of the numerical 
results on the interfaces until some of the analytic arguments are available for comparison. 



3.5.2 Asymptotic Analysis of the Interface Solutions: preliminaries 

Before addressing the large-w and small-u limits separately, let us put the TDGL equations into a form 
convenient for analysis and derive expressions for the length scales deep in the superconducting and normal 
regions. The disparity of these length scales in the large-u limit will motivate the boundary-layer analysis in 
that regime; while an inequality they satisfy will lead to the conclusion that J* — > Jc in the small-u limit. 
We make the substitution ■0 = /e*^, which yields 



uft = - f{9.)' + f-f, 

ui0t + ^i)f = 2fj, + fe^,, 

J = f'^Ox — Ma- 



(3.28) 
(3.29) 
(3.30) 



Next we restrict our attention to stationary solutions. Note that only spatial derivatives of appear now, 
allowing us to work with the superfluid velocity q — 9x instead of 6. The equations become 



fxx -q^f + f - f 
J* = fq-f^x, 



0, 



(3.31) 
(3.32) 
(3.33) 



where J* replaces J as these equations apply to the stall situation. Next multiply Eq. ( 3.32 ) by / an d note 
that the right hand side is now {f'^q)x which we can express in terms of /i by differentiating Eq. (3.33); these 
steps lead to 



Now let us assume the following asymptotic forms as a; — s- — oo: 

lim J{x) = /oo-/ie^/^M 



lim q(x) 

C — * — CO 

lim /i(x) 



qi e 



x/\g 



x/\fi 



(3.34) 

(3.35) 
(3.36) 
(3.37) 



Substituting these expressions into Eqs. ( 3.31 ), ( 3.33 ) and ( 3.34 ) and recalling the boundary conditions 
yields (is label "approach" correct?) 



^ ^le-^A. = 0, 



0. 



(3.38) 

(3.39) 
(3.40) 



Eq. ( 3.40 ) provides an expression for A^, the electric-field screening length. Since foo is always of 0(1), we 



see that shrinks as u — + 00 and diverges as w — > 0, which is consistent with the behavior seen in Figs. 3.8 



and 3.E 



More than one decay length appears in Eqs. (3.38) and (3.39). If they are not equal, the term with the 



Landau- Ginzburg Systems February 1, 2008, 10:21 



56 



shorter length is exponentially small compared to the other (s) and will not contribute to the x — > — cx) limit. 



Since none of the terms in Eq. (3.39) can equal zero individually, we conclude that the longer two of A/, Xq 



and A^ must be equal. Next, because the term multiplying e^^^' in Eq. ( [3.38| ) cannot equal zero on its own. 



we determine that Xq < Xf, making A^^ one of the longer lengths. Finally, if we assume that Xj 



Xu > Xa 



we find that A/ = 2'^/'^/'^ and A^ 



-1/2 



and reach a contradiction (except at u = 2). Thus provided 



the original assumption of an exponential approach is valid, we conclude that 

Af = A„ > A,,. 



(3.41) 



This equality of A/ and Xq is reasonable given that both / and q are related to the complex order parameter 
tjj. Also, having A/ > A^ is consistent with the large-w data seen in Fig. 3.8. If A^ 7^ A/ then 



A^ 



6/0 



A 



LA- 



(3.42) 



We identify this length scale as Xla since it coincides with that occurring in the solution of Eqs. ( ^.5.2| ) 
without any electric field (/^(a;) ~ 0), 



/'(a;) = /^-(3/^-2) sech^ 




(3.43) 



which was found by Langer and Ambegaokar in their study of phase slippage. The asymptotic form of 
Eq. ( 3.43| ) looks like Eq. ( 3.35| ) with A/ given by Eq. (3.42). As a matter of fact because A/ ^ A^ in the 
large-u limit, the profile of f{x) is only imperceptibly different from the Langer- Ambegaokar (LA) solution 
in the superconducting region and deviates from it only in the boundary layer, as is shown in Fig. |3^ . 

Recall that A^ diverges as m ^ 0; the inequality Ay > A^ implies that Xf must diverge as fast or faster in 
this limit. This scenario is consistent with the small-u data shown in Fig. 3.9 in which f{x) and E(x) vary 
on long length scales. Eq. (3.42) suggests that a diverging A/ implies that fao — + -\/2/3 and in turn that 
J ^ Jc as u ^ 0, 



Eq. (3.42) suggests that a diverging A/ imphes that fao - 
which is also consistent with what is found numerically. 
In the other asymptotic limit, deep in the normal regime, tp is very small and hence the nonlinear terms 
in Eqs. (3.3.1) can be dropped as was done for the bumps in the s mall- J limit. The result is a complex 
Airy equation, the asymptotic analysis of which was supplied in Eq. ( 3.21 ), where we saw the length scale 



Xav 



(uJ*) Somewhat like A^, XAiry shrinks as u 



00 and expands as u — > but with different 
iry, in the large- it limit motivates 



A^ and Xav 



powers of u. The presence of the disparate length scales, Xf 
the use of the boundary-layer analysis that comes next. We will see that XAiry scales in the same way as the 
boundary-layer thickness. 



3.5.3 Asymptotic Behavior of the Stall Current as ti — > cxd 



We have already seen in Fig. 3.8 that the large-w profile can be divided into two regions — one slowly varying, 
one rapidly varying, also known as the outer and inner regions, respectively. Furthermore, it has been 
suggested that the ratio of the length scales characterizing these regions decreases as u — > 00. These features 
make the problem ideally suited for boundary-layer analysis, in which one identifies the terms that dominate 
the differential equation in each region, analyzes the reduced equations consisting of dominant terms and 

then matches the behavior in some intermediate region. 

We start by eliminating the superfluid velocity q from Eqs. ( 3.5.2| ), resulting in 



/xx - {J* + ^ixfr'' + /-/' = 0, (3.44) 
f-xx - up II = 0. (3.45) 



Let us consider first the slowly varying, superconducting region. We saw in the preliminary analysis that for 
large u, /i(a;) is exponentially small, so we drop it. Next, let us assume that J* is small and drop it; we can 
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verify in the end that this is self-consistent. The reduced equation is 

f.. + f-f^ 0, 



(3.46) 



with solution f{x) = — tanh(a;/\/2). 

Moving in from the left toward the interface (into the boundary-layer region), / becomes small, and the 
second term in Eq. (3.44) which was subdominant becomes dominant. In this inner region / is small but 
rapidly varying, thus the dominant terms are 



(J* 



P 



(3.47) 



along with Eq. (3.45). Having identified the dominant terms, now we must make certain they balance. We 
assume that in the boundary layer, all the quantities scale as powers of u: 



f 



J* 



(3.48) 



Balancing terms in Eq. (3.47), we find 2a = j + S, while balancing terms in Eq. (3.45) yields 2{a + S) = 
1. Next, we need to ensure that the solutions in the boundary layer match onto the solutions in the 
superconducting and normal regions. By expanding the superconducting solution near the interface, we see 
that f{x) ^ —x/V2 as the boundary layer is approached; matching to the boundary layer requires fx ~ 1, 
so that a = 7. In the normal region, /x « —J*x, so that matching to the boundary layer requires fJ-x ^ J* , 
and /3 = J + S. Solving this set of equations, we conclude that a = j = S = 1/4 and /3 = 1/2, i.e., the 
stall current J* ~ u~^^^ for large u. Note that J* — > as w ^ 00, so that we were justified in dropping 

J- 1/4 j^^Q A^iry 



gives Xav 



-1/4 



, indicating that it may be 



J*V/3 from Eq. (^). Substituting J* - 
identified as the boundary-layer thickness. 

In order to determine the coefficient of the u~^^'^ in the stall current we need to solve the boundary layer 
problem. Let us rescale in the way suggested above: 



X = u-^/^X. 



(3.49) 



Substituting these rescaled variables into Eqs. ( 3.44 ) and ( 3.45 ), and then expanding F, M and J in powers 
of it~^/^, we obtain at the lowest order 



F, 



o,xx 



{Jo + Mo,x? 



0, 



Mo,xx - F^Mo = 0, 



with the boundary conditions (from the outer regions) 



i^o.x(-oo) =-l/V2, Mo(-oo) =0, 
Fo(-Hoo) =0, Mo,x(+oo) =-Jo. 



(3.50) 
(3.51) 



(3.52) 
(3.53) 



(As before we need an extra boundary condition to fix the translational invariance.) For an arbitrary Jq the 
solutions of Eqs. (3.5C) and (3.51) cannot satisfy the boundary conditions; Jq must be tuned to a particular 
value before all of the boundary conditions are satisfied, leading to a nonlinear eigenvalue problem for Jq. 
We have solved this eigenvalue problem numerically and find that Jq = 0.584491. Therefore, to leading order 
we have for the stall current 



J* « 0.584491 u"^/''. 



(3.54) 
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u 


J* 




V 




1 





3838 


0.3838 


0.01871 


0.01871 


5 





3407 


0.5094 


0.6400 


0.1914 


10 





3013 


0.5359 


1.573 


0.2797 


50 





2127 


0.5655 


8.258 


0.4315 


100 





1807 


0.5715 


15.59 


0.4931 


500 





1224 


0.5788 


62.51 


0.5875 


1000 





1033 


0.5807 


111.3 


0.6259 


5000 





0693 


0.5828 


407.9 


0.6847 


10000 





0583 


0.5833 


708.4 


0.7084 


50000 





0391 


0.5840 


2487 


0.7440 



Table 3.3: Representative numerical results for the stall ciurrent J* and kinetic coefficient r]. 
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This prediction agrees well with the numerical results and disagrees with Likharev's conjecture of a u~^^^ 
dependence |Q, as can be seen in Fig. 3.10| and in Table 3.3. It is in principle possible to carry out this 
procedure to successively higher orders, but the equations become cumbersome. Instead we have simply 
opted to fit our numerical data to a form inspired by the asymptotic analysis, 



J* = 0.584491 u"^/^- 0.117461 u"^/'^- 0.12498 u"^/'^ 
+0.163043 M"^/* + 0( u-"^!^ 



(3.55) 



Even with the higher order terms, the asymptotics are appropriate only for physically large values of u. 
3.5.4 Asymptotic Behavior of the Stall Current as ^ 

Now let us examine the opposite limit of m ^ 0. In this case the electric-field screening length becomes 
long, and Ivlev et al. [ |60[ have proposed that this makes the small-w limit useful for modeling gapped 
superconductors. As already suggested the inequality of length scales. Ay > implies that J* Jc- We 
will begin our small-w analysis by putting this result on firmer ground and extracting as a byproduct the 
u Q limit of the interface profile. 

The resettled equations. Recall that deep in the superconducting region A^ ^ u^^/"^. This observation 
suggests that we rescale distance: x — u^^/'^x; furthermore, to ensure that the normal current (— /ix) scales 
in the same way as the total current we rescale /i — u~^/'^jl as well. These rescalings yield 



ufxx - q" 

fi'f = "^qfx + fqx, 
J* = pq- 



0, 



(3.56) 
(3.57) 
(3.58) 



placing the small parameter u in front of fxx ■ If we expand these functions as series in powers of u 



f 

q 



= Mo 



/o + ufi + 
qa + uqi + . 
ufii + 



J* = J, 







U Ji* 



(3.59) 
(3.60) 
(3.61) 
(3.62) 



then we find at the lowest order 



-9o/o + /o - - 0, 
Mo/o = 2(7o/o,2; + foqo.x, 



(3.63) 
(3.64) 
(3.65) 



The solution of Eq. (3.63) is either fo — O (the normal phase) or /o = (1 — 9o)^^^ (the superconducting 



phase). Let us focus on the superconducting solutions. By eliminating qq, we obtain the first order equations 



lO.x 



Because /o ranges from foo to and fa 



/oVl - fo Ao 



Ao,s — /{ 



'2/3, we know that /o either starts at or passes through 




(Strictly speaking we should be writing here /oo,0j the lowest order term in the expansion of foo-) Thus, 



the effect of the pole in Eq. (3.66) must be considered. If it is not canceled by a zero in /io, fo,x diverges at 
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fo = 72/3. 

We can obtain an expression for /io(/o) by dividing Eq. ( 3.67 ) by Eq. ( 3.66 ), which leads to 

(2-3/0^) 



/to dfio = 



-dfo- 



(3.68) 



Integrating both sides and recaUing the boundary condition ^oo = 0, we find 

3 



2 ~ JO ^ Jo Joo 



+2 J* In 



-2 Jq* In 



4 •'"^ 



fo 



3^Vl-/o' 



l + y/T^ 

foo 



+ 3j*yw; 



2 

oo ' 



(3.69) 



where Jq = f^\/l^^f^- To keep /o,x from diverging, we insist that /io(/o = ^^2/3) = which can be shown 

'2/3, i.e. the small-w limit of the stall current is the critical depairing current. 

-co). 



from Eq. (3.69) to imply /o, 



Note that the pole in Eq. (p.66[ ) and the compensating zero in /io(/o) occur at the boundary {x 
We can rearrange Eq. ( ^.66| ) as follows 

(2 ~ 3/^) df _ . 
/o(o) /v^l - P Ao(/) 



(3.70) 



Then we can substitute in Eq. ( 3.69| ) for fLo{f), numerically integrate the resulting expression and finally 
invert it in order to calculate fo{x), the u — > profile. Figure |3.9| includes a comparison of /o(i) and the 
profile of a small-u numerical solution. 

To find the asymptotic behavior of fo and jlo in the superconducting region, Taylor expand /io(/o) around 

/oo 



Ao(/o) = -3V2(/o- 



(3.71) 



Notice that /io(/o) is a second order zero, so that /o.x = 0, as it should at the boundary. As a consequence, the 
integral supplying the inverse profile, Eq. ( 3.70 ), has a pole; integrating the expression in its neighborhood 
yields \/61n(-\/2/3 — fo), leading to 



foix) ~ ^2/3 - Ao exp (^x/Vtj 



(3.72) 



where Aq is an integration constant undetermined because of the translational invariance. Note /o(i) has 
the form assumed in the preliminary analysis with A/^o = \/u/6. Putting this result into Eq. (3.66) leads to 



fioix) 3V2AI exp (2x/V6^ 



(3.73) 



where A^^q — u ' /oo in agreement with the expression found previously. 

Let us examine Eqs. (3.66) and ( |3.67 ), which are strictly speaking superconducting solutions, in the 
normal (small- /o) limit. Eq. (3.67) leads to /io(^) ~ —JcX, and inserting this into Eq. (3.66) reveals that 
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/o — > in the following way 



fo{x) - exp {-JcX^/£) 



(3.74) 



This same dependence was seen earlier in the analysis of the bump shapes in the small- J limit, Eq. ( p.2C| ). 

What is surprising here is that what are ostensibly the "outer" equations for the superconducting region 
also satisfy the boundary conditions in the normal region and interpolate in between. This is consistent with 
the numerical observation that there docs not seem to be a boundary layer at small u, that the uf^x term is 
apparently not a singular perturbation. With this in mind, we pursue the perturbative expansion to higher 
orders. 

The 0{u) equations. The eigenvalue Jq was determined by examining the behavior deep in the super- 
conducting region and did not require imposing the boundary conditions on the no rmal side. Furthermore, 
the spatial dependence of the solution in this region is of the form assumed in Eqs. ( 3.5.2 ). We exploit these 
features to obtain higher order terms. The 0{u) equations are 



fo,. 



2go/o<7i - foh - 3/o'/i = 0, 



Ao/i + /oAi = 2(7o/i,£ + 2/04(71 + foqi.x + m,xfi 

Ji = 2/ogo/i + foQi 



The asymptotic form of fo{x) is 



(3.75) 
(3.76) 
(3.77) 



(3.78) 



and similarly for qQ{x) and fiQ{x). Eqs. ( 3.5.4 ) can be satisfied if the asymptotic form of fi{x) is 



(3.79) 



and similarly for qi{x) and fli{x). At 0{u^), f2(x) would have second-order polynomials multiplying the 
exponentials, and so on. Substituting these expressions into the differential equations allows us to determine 
the unknown constants (except for those associated with the translational invariance). For f^o it yields the 
series 



/I M 
^°°"V3 + 24^6 "^^768^6 



(3.80) 



(3.81) 



which corresponds to 

J* = 

Note that the first correction to the u 

The series found through the asymptotic perturbative expansion above can be obtained by another 
method. Looking back at Eqs. (3.72) and (3.73), we note that the ratio of decay lengths A//Ap = 2. If we 
insert the expressions we have for these length scales, Eqs. (3.4C) and (3.42), we find as u — > 



3V3 576%/3 5184V3 
limit of J* is of O(m^), since the lowest term is at the maximum 



hi 



-,1/2 



M 



2 „4 



(3.82) 
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0.010 




Solving for /oo , and then calculating J* , we find 

J* = J,(l - u/8)i/2(i _ m/24)-3/2, 



(3.83) 



with Jc — ■\/4/27, which when expanded for small-u agrees with the series (3.81) found above. We plot the 
small-M numerical data and this expression together in Fig. 3.11, The fit is surprisingly good at small u, 
suggesting to us that the corrections to Eq. ( |3.83| ) are exponentially small as m — > 0. 



3.6 Moving Interfaces 

At currents other than J* , the NS interfac es mo ve with a constant velocity. For such solutions the operator 
dt can be replaced by —cdx, so that Eqs. ( 3.5.2|) become 



u{-cq + fi)f = 2/j,g + fq^, 
J = fq- fJ-x- 



(3.84) 
(3.85) 
(3.86) 



While the boundary conditions on / and q remain the same, that on the scalar potential becomes ^oo = cgoo. 
Actually, it is more convenient to use instead p. = fi — cq, which is the gauge-invariant potential in the 
constant- velocity case. 

The superconducting phase invades the normal phase if J < J* and vice versa if J > J*. For currents 
near J*, the interface speed is proportional to (J — J*). In this linear response regime, one can define a 
kinetic coefficient (which Likharev refers to as a "viscosity" ) 



dc 
dJ 



(3.87) 



Figure 3.12 shows the numerically determined kinetic coefhcient as a function of u. For large-u, we find 
77 ~ for which we provide an argument below. 
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10' 10° 10^ 10^ 10^ 10-^ 10* 10^ lO'' 



U 

Figure 3.12: A log-log plot of the numerically determined kinetic coefficient as a function of u (the 
solid line) along with an asymptotic fit of 0.797u^^'^ (the dotted line). 



The velocity of the interface does not depend on the direction of the apphed current. This can be seen 
from the symmetry of the current in the equations we use, but it seems not to make physical sense. In fact an 
article by Gurevich and Mints [|4| addresses effects of the collision of quasiparticles with the NS boundary, 
and they do cause a small asymmetry in the normal zone propagation. The TDGL do not, however, account 
for kinetics across a phase boundary, so the effect is not included in Likharev's equations. 

Farther from the stall current, the velocities deviate from this linear behavior, as seen in Fig. 3.13. The 
greatest departure occurs in the limits J —^ and J — > Jc- In fact, Likharev |Q conjectured that the 
interface speed diverges in both of these limits; we find that it is bounded. 

The J — !■ limit. The moving interface equations, Eqs. simplify in the J ^ limit, since that 

limit implies that both g — > and /i — > 0, leaving only 



+ ucf^ + ./-/' = 0. 



(3.88) 



If we replace uc in the above equation by a speed v, then we have the steady-state version of Fisher's equation 
p9| , which is known to have propagating front solutions with w = 2 1^. In our case this implies that as 
J ^ 0, c = 2/m, which is in good agreement with the numerical data shown in Fig. 3.13 . 

We can combine the above result with an earlier one to suggest that 77 ~ u'^/* as u 00. In the large-u 



limit, we have information on the following two points: (1) the stalled interface (J — J* 
and (2) the interface in the absence of current ( J = 0, c = 2/u). In going from (1) to (2) 



-1/4 



current and velocity are A J 
might be approximated by 



'1/4 and Ac 



0); 

the changes in 
As M ^ 00, both of these changes are small so that ?/ 



77 



AJ 
A^ 



,3/4 



(3.89) 



yielding the behavior seen in the numerical data (see Fig. 3.12 and Table 3.3). 

The J —I- Jc limit. The numerical work indicates that the velocity is finite as J 
is shown in Fig, 



Jc, the limiting velocity 

3.14| as a function u. We can find an analytic bound on this velocity as follows. First, take 

Then 



Eqs. (p.6|), use the gauge-invariant potential fl, and find the constant- velocity analog of Eq. (B.34) 



Landau- Ginzburg Systems February 1, 2008, 10:21 



64 




substitute the asymptotic forms, Eqs. (3.5.2), into the resulting equations, leading to 



(^/^-A;2)^le-/^^=c,lA-VA^ 



2/oo9oo/ie 



x/\f 



if 



-MiA^^e^/^" = 0. 



(3.90) 
(3.91) 

(3.92) 



Arguments similar to those following Eqs. (3.38-3.40) lead one to the conclusion that in this case A/ = A^ = 



Ap. The above equations can then be shown to yield the following relation 



2 2 

u c 



+ [2 {uflX^ - 1) (3/^ - 2) - ufl + A"2] = 0, 



(3.93) 



where we have used = 1 — We find the bound by: (1) solving Eq. ( |3.93| ) for c; (2) substituting in 
foa = \/2/3 (which corresponds to J = Jc); and (3) extremizing that result with respect to the decay length 
A. The small-M limit of the resulting bound is — V2m/9, and the large-w limit is — 1/2-\/3. The square-root 
dependence of the velocity in the small-u limit agrees with the data. Now we can consider going from the 
stall current (J*,c = 0) to the critical depairing current (Jc,c ^ 
and Ac ~ u^/^, suggesting that the small-u kinetic coefficient 77 
the numerical data. We have also observed that as a function of J the speed appears to approach its bound 
via a square root dependence 



v}^'^) which results in changes AJ ~ 
' u'^/^, which is in rough agreement with 



c( J) = A + B{Jc - J) 



1/2 



(3.94) 



for all u. 
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3.7 Summary and Remarks 

In this chapter we have studied in detail the nucleation and growth of the superconducting phase in the 
presence of a current. The finite amphtude critical nuclei grow as the current is increased, with the amplitude 
eventually saturating as the stall current J* is approached, leading to the formation of interfaces separating 
the normal and superconducting phases. The stall current can be calculated in the limit of large u using 
matched asymptotic expansions, demonstrating once again the utility of this technique for problems in 
inhomogeneous superconductivity. We have also derived an analytic expression for the stall current for small 
u, which we believe to be correct up to exponentially small corrections. Deviations from the stall current 
cause the interfaces to move, and we have calculated the mobility of these moving interfaces for a wide range 
of u. Finally we have shown that the interface velocity c = 2/u as J ^ and that c is bounded as J Jc, 
in contrast to some conjectures in the literature. 

As in the magnetic-field analogy, the issue of stability and dynamics of the current-induced NS interfaces 
will be more complicated and interesting in the two-dimensional case. Some preliminary work in this direction 
has been reported by Aranson et a/.|^], who find that the current has a stabilizing effect on the NS interface. 
This can be interpreted as a positive surface tension for the interface, due entirely to nonequilibrium effects. 
They provide a heuristic derivation of an interesting free-boundary problem for the interfacial dynamics (a 
variant of the Laplacian growth problem); however, this free-boundary problem is sufficiently complicated 
that they were unable to solve it to compare with their numerical results. Clearly, further work in this 
direction would be helpful in understanding the nucleation and growth of the superconducting phase in 
two-dimensional superconducting films. 



Appendix A 

Detailed Derivation of the TDGL 



A.l Derive TDGL from Landau- Ginzburg Free Energy 

This note shows the details of deriving the time-dependent GL equations from the GL free energy. 

d d 



dt 5i)* 



d 

dip* 



Start with the GL free energy 



Now work out the functional derivative with respect to V'* 



G 



e.,A 



d 
dip* 

_d d_ 



{ihV + -A)iP ■ {-inV + -A)il}*\ = -A • {ihV + -A)V 

c c J c c 

{ihv + ^A)v • {-inv + ^A)v* = -inv ■ {ihV + ^A)V 



Put it together to find 

6tp* 



2m \ c 

1 / c \ ^ 

lib + Plthl'^ip -\ iihV + -A] lb 

2m V c / 



(A.1) 

(A.2) 

(A.3) 
(A.4) 
(A.5) 
(A.6) 

(A.7) 
(A.8) 
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Now do the next independent variable, A. First we write out some simple definitions for reference. 



V X A : 
V X V X A : 



OA, dA 



dy dz 
d fdAy dA, 



dA, dA, 



dx 



dz 



J + 



dAy _ dA^ 
dx dy 





(dA,, 


dAA] 


dz 


[ dx 


dz J 



d fdAy dA, 



dx \ dx 



dy 





(dA, 


dAy\-\ 


) dz 




dz ) 



dx \ dx 



dA^A 


-A 


(dA, 


dAy\-\ 


dz ) 


dy 


\ dy 


dz ) 



(A.9) 
(A.10) 
(All) 
(A.12) 



Now lot's look at the functional derivative of (V x A)^. When we look at a single component, it becomes a 
little clearer. 



= 



d d{V X A)^ 
ax d^ 

OX 

d diV X Af d . , , J, 

dy d^ dy 

d d{V X A)2 d 



dz d'-t 

You can see we are reconstructing a cross product 



= -2(VxA).i 



(5(V X A)^ 



d{V X A)^ d{V X A), 



dy 



dz 



to get 



Similarly, 



6A 



(V X A)2 = 2V X V X A 



5A 



(V X A) • H = V X H 



Now figure out the rest of the derivative for A. 



_d_ 
dA 



{ihV + -A)ip ■ {-inV + -A)V'* = -tpi-ihV + -Al^)* + c.c. 

c c J c c 



We can add the two pieces together to find the current 

- — (V'VV*-V*V^) + 
2mc 

Putting the second equation together, we find 
5G 1 



_ = — V X (V X A - H) - ^ 

5A A-K ^ ' 2mc 



(A.13) 

(A.14) 
(A15) 

(A.16) 
(A.17) 
(A.18) 
(A.19) 
(A.20) 
(A.21) 
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A. 2 Find Dimensionless Variables 

We start now with the known form of the TDGL, 

^ + f c/-^) + aV' + /3|^l'^ + ^ (i W + -^Afi; = (A.22) 

+ aV(j) + —V X V X A - —N,Vi)* - f*Vij) + — A|t/.|2 = (A.23) 
c at 4tt 2m mc 

We begin by requiring ^jj vary between zero and one so that ip = ipoip' where ipo = y/—a/f3. That determines 
the overall factor on the first equation. 

We can immediately determine A = AqA' to be 



= \ ■ (A.25) 



We also see that the gradient's factor is of the form 

n 1 _ 1 _ ^ 
y/-2ma X K X 

where x = Ax'. 



(A.26) 



^ + ^hA +^'- IV'I V - (- V + A'fi,' = (A.27) 



—aV \ dt h ' J K 
Looking now to the second equation, we can divide through by Aq to see 



' ■ ' vV + -^V X V' X A' + 7r-n5\/ 7^ ^(V' vV - V' vV) ^A'|vf = o 



cto dt' V 2mac2 A AirX"^ 2mXP V 2TOac2 mc/3 

(A.28) 

The last term is the only one we already know. Multiplying all terms by the inverse of its prefactor gives 



^^TT + ^ I lo VV - - — TT^V X V' X A' - — \ - iih'V'ib * - V VV) + A' V r = 

cto e^a dt' K^elV'oP Aire^aX^ 2X]j 2ma^^ ^ ^ j ^ ir i 

(A.29) 

Now we know from the third term that we must define 



With that, the third term simplifies dramatically to yield 

+ £SF''*' + V X V X A' - ±»'V'i,'- - VV) + A'|,f = (A.31) 



Landau-Ginzburg Systems February 1, 2008, 10:21 69 
The last definitions look clear. 

to = — ^ and ,^o = Kh/'o A.32 

am 

If we return to the equation for i/» to finish changing its variables, we find 

a (dip' ie \ 1 ( d%l>' ie4>oto , \ r \ oo\ 



We are hoping to find that 

iecpoto 



= iK = i-, (A.34) 



n 

which it is. Our final equations, when this becomes clearer, will be (dropping primes) 

7 + iK<P^^ IV'I'^ + V + Af^ = (A.35) 

— +V?!> + VxVxA - — (V-VV* - V* V^) + A|^|2 = (A.36) 



Appendix B 



Likharev's Equations as an Active 
Kinetic Equation 



Likharev's equations behave like nonlinear diffusion equations. They don't appear to be the same at first 
glance, however. We would like to see whether they could resemble more traditional diffusion equations, so 
we begin with Likharev's equations 

uft - f..-fel + f-f (B.i) 
u{et + ^i)f = ife,), (b.2) 

J = fOx-fix- (B.3) 



and try to coax them into a form like that discussed in Gurevich and Mints |5q ] 

= ^S-m^,/?) (B.4) 



where /3 is a parameter representing our {u,J). Wc may be able to put our equations in this form if we 
change variables from (/, 0) to {f,js) where js = f^Ox is the supercurrent. Equation B.I is already in the 
correct form if we substitute 6^ — j / P ■ 

uft = h. + f{l-f~ ^) (B.6) 

If we take the derivative of |B.2|, we get 



uiOxt + M.)/' + "^jifOx) = {fex)xx (B.7) 
We can get rid of 9 by constructing the derivative of the supercurrent 

ifOx)t = 2fft0x + fO^t. (B.8) 

We substitute ft from |B.l| 

2ufftex = 2fdx{uft) = 2fdx{fxx -fel + f- /3) (B.9) 
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into the derivative of the supercurrent to find 

n{fO^)t = if 9^)^^ - 2^(/20,), + u{J - fe.,)f + 2f0,U, - 2fel + 2(1 - f)fe.,. (B.IO) 



It is time to substitute pO^ = j to find 

ujt = J.x " |(/xJx - /xxj) - (2 + + 2j + ujf - (B.U) 

We could condense this a httle to 



,^(^)^, 



ujt=U.-2'-[^] +u{J -j)p + 2j{l-f-^]. (B.12) 



What remains is of the form of an active kinetic system 



^=F.(xi,X2,...,x„) + |:(x:A,^| (^ = l,2,...,n) (B.13) 



as described in Vasil'ev, Romanovskii, and Yakhno |103|, but the numerous spatial derivatives complicate 



the first integrals typically used to examine such equations. 



Appendix C 

Amplitude of the Critical Nuclei in 
the J ^ Limit 



In this appendix we provide a self-consistent calculation of the amplitude of the critical nuclei in the J ^ 



limit. Choosing the gauge appropriate for bumps centered at x = and combining Eqs. (3.3.1) into one 
equation yields 



[-udt + iuJx + 9^ + 1] il;{x, t) = |?/;(a;, t)\^ %p{x, t) 



-lU 



dy Im {ilj*{y,t)dyilj{y,t)) 



ip{x,t). 



(C.l) 



The propagator for the linear operator appearing on the left hand side of Eq. ( p.l[ ) satisfies the condition 

(C.2) 



[-udt + iuJx + dl + \\ G{x, x; t - t')e(t - t') 

-u5{x-x')5{t-t') 



and is given by 



G(x, x' ] t) 



{—) 
\4tttJ 



1/2 



exp 



X exp 



iJT{x 



T 

u 12u ^ 
x') u{x — x'Y 



At 



(C.3) 



Ivlev et al. |9[ used this linear propagator to evolve perturbations having widths of 0(1) and carrying 
no current. Without the nonlinear terms such perturbations initially grow but ultimately reach a maximum 
size and then decay away. Ivlev et al. suggested that the amplitudes of the critical nuclei are exponentially 
small in the J — > limit by asking what sized initial perturbations are of 0(1) at their maxima. Their 
arguments motivated us to use the propagator in a more careful estimate of the amplitude that includes the 
nonlinear terms as an essential ingredient. 

We can convert Eq. (C.l) into an integral equation by multiplying both sides of Eq. (C.l) (with x — > x' and 
t t') by G{x, x';t — t') and integrating over all x' and integrating t' from to t. After some manipulation 
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these steps lead to 



dt' 



dx' G{x,x',t~t') 



iPix', t')S{t - t') - -\i,{x' ,t')\^ ^(x', t') 



dy im{riy,t')dyi'iy,t')) 



^{x',t')}, 



(C.4) 



where t > 0. 

In order to estimate the amphtudc of the threshold solutions, we will substitute into Eq. ( p.4[ ) the 
following form 



ipix,t) = ?Aoexp 



uJx^ 



(C.5) 



Note that this fo rm is stationary and has a fixed Gaussian shape (which is inspired by our WKB approxi- 
mation, see Eq. ( ^.20| )) but it has an arbitrary amplitude which we will determine self-consistently. 

Let us take the t ^ oo limit and focus on x — Q since our interest is in the amplitude. After substituting 
Eq. ( |C.5| ) into Eq. ( |C.4| ), we can do both integrals for the first term on the right hand side exactly, and it 
can be seen to decay to zero in the t oo limit. Next, we perform the x' integration of the second term on 
the right hand side (//), which yields 



dr 



uJ Jo VI + 3' 



; exp 



24r2 - 4t3 - 3t'^ 
12uJ(l + 3r) 



where r = Jt. We now apply the method of steepest descent to obtain 



^/2^tJ 



exp 



32 



81 uJ 



(C.6) 



(C.7) 



In the third term on the right hand side (///) of Eq. (C.4), we make the substitution y = vx' and then 
perform the x' integration giving 



/// 



mJ2 



dr 



dv- 



(2r + t2 



X exp 



y/[l+T{l + 2v^)]^ 

24^2^2 _ _ _^ 2v^)t^ 
12uJ[l + (1 + 2t72)r] 



(C.8) 



The maximum of the term in the exponential of /// occurs at u = 1 (which is an endpoint). Linearizing 
about that maximum provides 



///; 



dTT{2 + t) r 24t2 - 4r3 - Sr" 



dw exp 



(C.9) 



uJ(1 + 3t)2 

where w = 1 — v. After the w integration, we apply the method of steepest descent to the r integration to 



Landau- Ginzburg Systems February 1, 2008, 10:21 



74 



obtain 



///: 



TTUIpQ 9 



2J 



■ exp 



Putting all of these results back into Eq. (C.4) gives 



■ exp 



32 



81 uJ 



f 32 1 


■9 1" 


[81 uJ ) 


8 ^ u_ 



(C.IO) 



(c.ii) 



which provides the expression given in the text, Eq. (3.22). This calculation clearly runs into trouble when 



u < 8/9; however, the numerical coefhcients in front of these integrals are sub-leading terms, and they can 
be varied by adding sub-leading terms to the initial Gaussian guess. 
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